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NOTATION 

r~\  The  horizontal  mean  (often  just 

capitalized:  g  =  T,  etc.) 

<  >  The  vertical  mean  over  the  mixed  layer 

b  =  B  +  b  Buoyancy 

B*  Parameter  for  relative  importance  of 

H*  to  Ro-1  for  cyclical  steady  state 

5^"  Turbulent  buoyancy  flux 

^7  Complex  turbulent  momentum  flux 

r  Specific  heat  at  constant  pressure 
CP 

C  Drag  coefficient  in  quadratic  wind 

D  stress  law 

C  =  U  +  iV  Mean  horizontal  velocity  in  complex  form 

c  =  u  +  iv  Turbulent  horizontal  velocity  in  complex 

form 

d  Ekman  depth  of  frictional  resistance 

d1  Surface  production  zone  (Niiler);  also 

used  as  a  depth  proportional  to  Ekman 
depth  (d) 

D  Total  rate  of  viscous  dissipation  of 

turbulent  kinetic  energy  in  the 
mixed  layer 

D  Molecular  diffusion  coefficient 

s 

1  F  =  L,  m  =  u2  + v2  +  w2    Turbulent  kinetic  energy  (per  unit  mass) 

2  2ri  i       2 

E*  <E>/(hAB) 

*     *     * 

•11   +  E22   +  E33 


eTi  +  Eoo  +  Eqo  Turbulent  kinetic  energy  nondimension- 


alized  with  the  friction  velocity 
squared,  u2. 


f  =  21^1  sin*  Coriolis  parameter 


F  Damping  force  for  inertial  motions 

g  Apparent  gravitational  acceleration 

G  Total  rate  of  mechanical  production 

of  turbulent  kinetic  energy  in  the 
mixed  layer 

h  Depth  of  mixed  layer 

h  Mixed  layer  depth  if  the  vertical 

r  extent  of  the  region  is  stationary 

or  retreating 

H  Mixed  layer  depth  (h)  non-dimension- 

al i  zed  with  the  Obukhov  length 

h  Compensation  depth  for  radiative  heating 

J  Mechanical  equivalent  of  heat, 

4.186  x  107  9m  cm2 
cal  sec2 

K0  Eddy  mixing  coefficient  for  neutral 

conditions 

K  =  K  Eddy  viscosity 

K  Eddy  conductivity 

L  Obukhov  length  scale 

Integral  turbulence  length  scale 

o  Convective  length  scale 

i  Rotational  length  scale 

m.  ,  n.  Model  constants 

i  '  Ki 

n  Value  for  integrated  flux  Richardson 

number 

N  Brunt-Vaisala  frequency 

p  Fluctuating  pressure  component 

gp  /ax  Geostrophic  pressure  gradient  in  the 


g   a 


ocean 


VI 


PE  +  KE 

* 
P 

Q 

Q' 

Qo 


Total  mechanical  energy  per  unit  mass 
in  the  mixed  layer 

Nondimensional  entrainment  flux 

Radiation  absorption 

Solar  heating  function 

Solar  radiation  through  interface 
(incident  minus  reflected) 


1  n2   u2  +  v2 

2  q  ="~ 2 


Ro 
Ri 


RF 


Horizontal  component  of  turbulent 
kinetic  energy 

Fraction  of  net  production  going  to 
potential  energy  by  means  of 
entrainment 

Rossby  number 

Gradient  Richardson  number, 

9B/8Z 

|3C/3Z|2 

Flux  Richardson  number, 


-  bw/  (^  fy  +  vw  |^) 


Ri 


Ri 


Ri. 


s  =  S  +  s 
t 


Gradient  Richardson  number  in  the 
entrainment  zone 

Overall  or  bulk  Richardson  number, 
hAB/ | AC | 2 

Total  Gradient  Richardson  number 

9b/9z 
(au/8z)2 
Salinity 

Time 


VI  i 


(u.)  =  (U+u,  V+v,  W+w) 


U9'V9 


u*b* 


W 


(x..)  +  (x,y,z) 


Total  instantaneous  velocity,  the  sum 
of  the  mean  and  turbulent  components 
(neglecting  any  geostrophic  com- 
ponent) 

Geostrophic  components  of  velocity 

dh/dt,  entrainment  velocity 


Friction  velocity,  v|cw(o)| 

Surface  buoyancy  flux,  -bw(o) 

Wind  velocity;  sometimes  used  as  mean 
vertical  water  velocity 

Nondimensional  entrainment  velocity, 


u  /V<E> 
e 


Rectangular  space  coordinates  with 
x.  =  z  aligned  vertically  upward, 

parallel  to  the  local  apparent 
direction  of  gravity 

-Pq-^p/sS 
Coefficient  of  thermal  expansion, 


Pn_13p/3T 


(z) 


'T 
6' 
i 
5C0 

/.B 


Extinction  coefficient  for  net  solar 
radiation 

3B/3Z  for  Z  <  (-h-6)  =  N2 

9T/3Z  for  Z  <  (-h-6) 

Fraction  of  turbulent  energy  dissipated 

Thickness  of  entrainment  zone 

Excess  surface  velocity,  C(z=o)  -  <C> 

Change  in  mean  buoyancy,  across  the 
entrainment  zone,  B(-h)  -  B(-h-6) 


vi  n 


AC! 


Mean  velocity  drop  across  the  density 
interface 

Rate  of  viscous  dissipation  of  tur- 
bulent kinetic  energy, 


-  _   Qik(x  -  cnt) 

*  -  V    l 


9U.8U. 


Internal  wave  amplitude  on  interface 


=  T  + 


Temperature 


K 

A 

M 
v 

P 


p0=p(vv 


Temperature  flux  scale,  b  /eg 

Thermal  conductivity 

Heaviside  unit  step  function;  also 
used  as  dissipation  length  scale 

Molecular  viscosity 

Kinematic  viscosity 

The  instantaneous  density  of  the  sea 
water 

The  density  for  representative  values 
of  salinity,  s0,  and  temperature,  e0 


Representative  density  of  air 


p«C„ew 
o  p 


Turbulent  heat  flux 


1*1 

r 
e 


Surface  stress,  p„  u? 


Dissipation  time  scale 


Dimensionless  time  scale,  t/x 


TI 

4> 


Time  constant  for  damping  by  F 


Latitude 


IX 


*  c*/d-*    <E>  |aC|2 

E  /Ri   =  ' L— 

(hAB)2 


* 


H'  Nondimensionalized  turbulent  kinetic 

3 
energy,  h<E>  /  (u*  t  ) 

Earth's  rotation  vector 

a)  Angular  frequency 


A  GENERAL  MODEL  OF  THE  OCEAN  MIXED  LAYER 

USING  A  TWO-COMPONENT  TURBULENT  KINETIC  ENERGY 

BUDGET  WITH  MEAN  TURBULENT  FIELD  CLOSURE 


Roland  W.  Garwood,  Jr. 


A  non-stationary,  one-dimensional  bulk  model  of  a  mixed  layer 
bounded  by  a  free  surface  above  and  a  stable  nonturbulent  region 
below  is  derived.   The  vertical  and  horizontal  components  of  turbu- 
lent kinetic  energy  are  determined  implicitly,  along  with  layer 
depth,  mean  momentum,  and  mean  buoyancy.   Both  layer  growth  by  en- 
trapment and  layer  retreat  in  the  event  of  a  collapse  of  the  ver- 
tical motions  due  to  buoyant  damping  and  dissipation  are  predicted. 
Specific  features  of  the  turbulent  energy  budget  include  mean  tur- 
bulent field  modeling  of  the  dissipation  term,  the  energy  redistri- 
bution terms,  and  the  term  for  the  convergence  of  buoyancy  flux  at 
the  stable  interface  (making  possible  entrainment) .   An  entrainment 
hypothesis  dependent  upon  the  relative  distribution  of  turbulent 
energy  between  horizontal  and  vertical  components  permits  a  more 
general  application  of  the  model  and  presents  a  plausible  mechanism 
for  layer  retreat  with  increasing  stability.   A  limiting  dissipation 
time  scale  in  conjunction  with  this  entrainment  equation  results  in 
a  realistic  cyclical  steady-state  for  annual  evolution  of  the  upper- 
ocean  density  field.   Several  hypothetical  examples  are  solved,  and 
a  real  case  is  approximated  to  demonstrate  this  response.   Of  par- 
ticular significance  is  the  modulation  of  longer-period  trends  by 
the  diurnal-period  heating/cooling  cycle. 


1.   INTRODUCTION 

1.1  Purpose  of  the  Study 

The  objective  of  this  study  is  the  formulation  of  a  unified  mathematical 
model  of  the  one-dimensional,  nonstationary  oceanic  turbulent  boundary  layer. 
In  particular,  this  model  should  help  explain  and  predict  the  development  in 
time  of  the  seasonal  pycnocline. 

Interest  in  the  ocean  mixed  layer  stems  from  both  theoretical  and  prac- 
tical considerations.  Thermal  energy  and  mechanical  energy  received  from  the 
atmosphere  not  only  control  the  local  dynamics,  but  the  layer  itself  modu- 
lates the  flux  of  this  energy  to  the  deeper  water  masses.  Conversely,  flux 
of  heat  back  to  the  atmospheric  boundary  layer  has  an  important  influence 
upon  the  climate  and  its  fluctuations.  Figure  1  depicts  the  mechanical 
energy  budget  for  the  ocean  mixed  layer. 
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In  addition  to  climatological  and  ocean  circulation  studies,  other 
applications  of  a  practical  model  of  the  mixed  layer  include  investigations 
of  advection  and  turbulent  dilution  of  dissolved  or  suspended  concentrations 
of  matter  such  as  pollutants  and  nutrients.  Prediction  of  entrainment  of  the 
deeper,  nutrient-rich  water  into  the  mixed  layer  has  particular  pertinence  to 
the  estimation  of  primary  productivity. 

The  need  to  further  develop  and  attempt  to  improve  the  one-dimensional 
mixed-layer  model  is  evident.  The  failure  of  earlier  models  to  consistently 
explain  the  annual  cycle  of  thermocline  development  is,  for  the  most  part, 
not  attributable  to  the  assumption  of  one-dimensionality. 

A  refined  treatment  of  the  often  neglected  terms  of  the  turbulent  kin- 
etic energy  budget  promises  to  improve  the  physical  understanding  of  the 
turbulent  ocean  boundary  layer  and  make  possible  the  creation  of  a  better- 
performing  model . 

A  detailed  review  of  the  literature  associated  with  mixed-layer  modeling 
is  undertaken  in  Chapter  II,  but  Table  1  summarizes  the  major  historical 
contributions  to  the  understanding  of  the  physics  of  the  mixed  layer. 

The  works  of  primary  concern  here  are  those  dealing  explicitly  with 
equations  for  the  production,  alteration,  and  destruction  of  turbulent  kin- 
etic energy  within  the  mixed  layer.  Kraus  and  Turner  (1967)  were  the  first 
to  heed  the  turbulent  kinetic  energy  budget  in  a  one-dimensional  mixed  layer 
model,  utilizing  the  approximately  decoupled  state  of  the  equations  for  the 
thermal  and  mechanical  energies.  By  neglect  of  the  frictional  generation  of 
heat,  the  vertically  integrated  heat  equation  becomes  a  relationship  for  the 
conservation  of  potential  energy.  However,  viscous  dissipation  cannot  be 
neglected  in  the  turbulent  energy  budget,  and  Geisler  and  Kraus  (1969)  as 
well  as  Miropol'skiy  (1970)  and  Denman  (1973)  added  it  to  the  basic  model. 
Niiler  (1975)  showed  that  in  addition  to  the  equations  for  thermal  (poten- 
tial) and  turbulent  (kinetic)  energy,  an  equation  for  the  mean  kinetic  energy 
should  properly  be  incorporated  because  entrainment  converts  some  of  the  mean 
flow  energy  into  turbulent  energy,  over  and  above  the  parameterized  wind- 
stress  production. 

Further  questions  remain  and  limit  the  general  applications  of  these 
earlier  models: 

a.  The  viscous  dissipation  of  turbulent  energy  has  been  parameterized 
as  a  fixed  fraction  of  the  wind-stress  production  and  hence  is  a  function  of 
the  surface  friction  velocity  u*.  Dissipation  may  be  related  to  the  integral 
velocity  scale  of  the  turbulence,  but  this  scale  is  not  always  proportional 
to  u*.  Surface  heat  (buoyancy)  flux  and  entrainment  fluxes  can  contribute 
significantly  to  the  turbulent  intensity. 

b.  Entrainment  may  also  be  considered  a  function  of  the  ambient  turbu- 
lence parameters.  Instabilities  leading  to  entrainment  are   probably  induced 
by  horizontal  turbulent  velocities  locally  at  the  bottom  of  the  mixed  layer, 
so  the  entrainment  rate  doesn't  necessarily  correspond  to  an  integral  con- 
straint upon  the  total  turbulent  energy  budget. 
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Date 


Author(s) 


Cont r i but  ion 
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Anderson 
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Ki  ta  igorod 

1961 

Kraus  and 

Rooth 

1967 


1969 


970 


1973 


197^ 


Kraus  and 

Turner 


Ge  i  s  ler 
and  Kraus 


Mi  ropol ' skiy ; 
Denman 


Pol  lard  ,  Rhi  nes 
and  Thompson 


Ni  i 1 er 


Constant  eddy  viscosity  solution  to  steady-state  momentum 
equation:  Ekman  spi  ral ;  "depth  of  frictional  influence"; 
suggested  h  *  W//sin  tji  . 


Improved  current  measurements  to  demonstrate  h 
(Rossby  number  ~  constant). 


u,/f 


Simultaneous  solution  of  steady-state  heat  and  momentum  equa- 
tions using  eddy  mixing  coefficients  variable  with  Richardson 
number. 


Using  dimensional  analysis,  suggested  that  h 
length  scale. 


L,  the  Obukhov 


Penetration  of  solar  radiation  to  depth  makes  steady-state 
possible  for  heat  equation  because  of  surface  heat  loss  by 
conduction,  evaporation,  and  long  wave  radiation;  unstable 
density  profile  above  compensation  depth  (hc)  is  source  of 
turbulent  kinetic  energy  produced  by  convection. 

Included  mechanical  stirring  (parameterized  in  terms  of  wind 
stress)  as  well  as  convective  production  as  important  source 
of  turbulence  for  mixing:  turbulent  kinetic  energy  equation 
and  heat  equation  form  two-equation  model  in  two  unknowns--T, 
h;  non-steady-state  solutions  and  "retreating"  (h)  possible; 
viscous  dissipation  neglected. 

First  "slab"  model  in  which  momentum  equation  is  solved  to- 
gether with  turbulent  energy  equation  and  heat  equation;  layer 
assumed  homogeneous  in  T  and  U,  V  and  hence  moves  as  a  slab; 
assumed  buoyancy  flux  is  fixed  portion  of  mechanical  produc- 
tion—essentially same  as  Kraus  and  Turner;  applied  model  to 
atmosphere  with  subsidence. 

Assumed  dissipation  a  fixed  fraction  of  mechanical  production; 
remaining  turbulent  kinetic  energy  goes  to  buoyancy  flux  down- 
wards including  entrainment;  essentially  same  as  Kraus  and 
Turner. 

Slab  model  applied  to  ocean  mixed  layers;  used  total  mechan- 
ical energy  equation  (rather  than  turbulent  equation)  plus 
momentum  and  heat  equations;  h,  T,  U,  V--all  functions  of 
time.   Ignored  effects  of  turbulent  energy  budget  altogether. 

Re-emphasized  the  need  to  use  the  turbulent  kinetic  energy 
budget  apart  from  the  total  mechanical  energy  budget;  divided 
region  into  three  sub-regions  with  no  mechanical  production 
in  most  of  the  mixed  layer.   Included  turbulent  kinetic 
energy  produced  by  entrainment  of  zero-velocity  water  into 
moving  slab,  but  surface  mechanical  production  minus  dissi- 
pation still  parameterized  in  terms  of  u.,..   Dissipation  is 
not  affected  by  the  additional  entrainment  production,  pos- 
sibly causing  a  too-large  entrainment  rate. 


c.  The  use  of  the  total  turbulent  energy  equation  and  consequently  the 
neglect  of  energy  redistribution  among  components  also  results  in  a  somewhat 
inconsistent  method  of  predicting  layer  "retreat."  The  consideration  of 
separate  budgets  for  the  horizontal  and  vertical  turbulent  energy  components 
will  permit  a  more  consistent  interpretation  of  both  entraining  and  retreat- 
ing mixed  layers. 

In  this  paper,  ad  hoc  mean  turbulent  field  modeling  of  the  terms  in  the 
turbulent  kinetic  energy  equations  permits  the  inclusion  of  these  often 
salient  effects  in  a  generalized  one-dimensional  model  of  the  ocean  mixed 
layer. 


2.  CHARACTERISTICS  OF  THE  OCEAN  MIXED  LAYER 

The  ocean  mixed  layer  is  defined  as  that  fully-turbulent  region  of  the 
upper  ocean  bounded  on  top  by  the  sea-air  interface.  The  wind  and  intermit- 
tent upward  buoyancy  flux  attributable  to  surface  cooling  are  the  primary 
sources  of  mechanical  energy  for  the  mixing. 

The  most  distinctive  feature  of  this  layer  and  what  really  defines  its 
extent  is  its  relatively  high  intensity  of  continuous,  three-dimensional 
turbulent  motion.  Vertical  turbulent  fluxes  within  the  mixed  layer  can  be  so 
much  greater  than  vertical  fluxes  in  the  underlying  stable  water  column  that 
the  dynamics  of  the  layer  are  essentially  decoupled  from  the  underlying 
region.   (Of  course  the  dynamics  of  the  underlying  water  masses  are  probably 
\/ery   dependent  upon  the  mixed  layer.) 

Typically,  an  actively  entraining  mixed  layer  is  bounded  on  the  bottom 
by  a  sharp  density  discontinuity  separating  the  layer  from  a  stable,  essen- 
tially nonturbulent  thermocline.  Minimal  stress  at  the  bottom  together  with 
high  turbulence  intensity  results  in  an  approximate  vertical  uniformity  in 
mean  velocity  and  density.  This  ostensible  homogeneity  is  the  root  of  the 
term  "slab,"  often  used  to  describe  the  layer.  On  the  other  hand,  only  small 
gradients  in  these  mean  variables  give  rise  to  large  turbulent  fluxes.  There- 
fore, even  the  slight  non-homogeneity  is  highly  important  in  the  physics  of 
the  region  and  should  not  be  neglected  at  the  outset. 

The  nearly  zero-flux  state  of  the  underlying  thermocline  causes  the 
bottom  boundary  condition  of  the  mixed  layer  to  act  almost  as  a  slip  condi- 
tion on  the  mean  velocity.  This  in  turn  creates  a  trap  for  inertia!  motion. 

Deepening  of  the  mixed  layer  is  accomplished  by  entrainment  of  the  more- 
dense  underlying  water  into  the  turbulent  region  above.  This  process  entails 
a  potential  energy  increase  and  cannot  take  place  without  an  energy  source-- 
the  turbulent  kinetic  energy  of  the  mixed  layer  above. 

A  simplified  picture  of  the  region  is  shown  in  Figure  2.  There  is  an 
appealing  practical  aspect  to  the  judicious  use  of  the  assumption  of  vertical 
homogeneity.  This  assumption  permits  the  use  of  the  vertically  integrated 
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momentum  and  heat  (buoyancy)  equations,  thus  avoiding  the  turbulent  flux  (of 
momentum  and  buoyancy)  closure  problem  altogether. 

The  depth  of  the  oceanic  wind-mixed  surface  layer  is  typically  on  the 
order  of  ten  to  a  hundred  meters.  The  horizontal  scale  size  is  that  of  the 
radius  of  the  circle  of  inertia—seldom  larger  than  a  few  kilometers  in  tem- 
perate latitudes.  These  two  dominant  scale  sizes  are  usually  significantly 
smaller  than  the  horizontal  scale  sizes  of  the  driving  meteorological  distur- 
bances, water  mass  features,  and  distance  to  lateral  boundaries.  Therefore, 
the  approximation  of  local  horizontal  homogeneity  for  all  mean  variables  is 
usually  accurate  and  is  a  basic  assumption  of  this  work.  The  local  conse- 
quence of  some  lateral  inhomogeneity  can  be  parameterized  without  qualita- 
tively undermining  such  a  one-dimensional  model.  For  example,  a  divergence 
of  the  horizontal  current  field  results  in  a  nonzero  vertical  mean  velocity 
which  in  itself  can  be  assumed  locally  uniform  in  the  horizontal.  A  minor 
rate  of  loss  (compared  with  the  surface  flux)  of  mean  momentum  by  lateral 
and/or  downward  radiation  of  inertial  motion  may  be  parameterized  also  with- 
out compromising  the  dominant  processes. 

Substantial  barotropic  and  baroclinic  features  in  the  mean  fields  can 
be  linearly  superimposed.  The  mean  fields  of  concern  are  therefore  the  hori- 
zontally homogeneous  components  of  the  total  fields.  In  particular,  the 


momentum  equation  has  the  geostrophic  component  subtracted  out,  eliminating 
any  lateral  pressure  gradient. 

1.3  Fundamental  Principles  and  Equations 

The  underlying  principles  employed  in  studying  the  mixed  layer  are  the 
combined  conservations  of  mass,  momentum,  thermal  energy,  and  mechanical 
energy. 

The  conservation  of  momentum  and  the  condition  of  incompressibility  are 
reflected  by  the  Navier-Stokes  equations  of  motion,  invoking  the  Boussinesq 
approximation: 

2 

9ui  aui  3~  3  ui 

p0  3t-+  po  "j   377  +    poeijkfi/k  +  377"   (po  -  p)9  6i3  =  ^377377    •      (K1) 

J  JO 


3U, 

Because  frictional  generation  of  heat  is  negligible  compared  with  typi- 
cal magnitudes  for  the  divergence  of  heat  flux,  the  conservation  of  thermal 
energy  is  decoupled  from  the  mechanical  energy  budget,  and  the  first  law  of 
thermodynamics  for  an  incompressible  fluid  gives  the  heat  equation. 

ii  +  ~    1L_  =  JL-     3  9       .  _Q_  (i  7\ 

»t       Uj   3xd  -p0Cp3xj3Xj       p0Cp         •  U"J' 

The  conservation  of  salt  mass  is  of  the  same  form,  but  it  lacks  a  term 
analogous  to  the  radiation  absorption  term,  Q. 

3S  ,  ~   3S   _  n    3  s  /-,  „\ 

3t  +  uj  3xTs  DS33T3T7-  (1'4) 

A  simplified  but  sufficiently  accurate  equation  of  state  for  local 
and  reasonably  short-term  application  in  the  mixed  layer  is  given  by 

p  =  Pn  [l  -  eCe  -  60)  +  a(s  -  s0)j  (1.5) 


where  p0  =  o'(e0,s0)  is  a  representative  but  arbitrary  density  at  the  time 
and  location  of  consideration.  The  coefficients  (b)  and  (a)  are  assumed  to 
remain  constant. 

Taking  the  scalar  product  of  (uf-f )  with  the  respective  terms  of  equation 
(1.1),  the  mechanical  energy  equation  is  formed, 


U  •  U  •    \  „        /U.U.\        U  •        „~ 

1  1  J   -   _3_  [   1  1  1.1   9£ 

2  /   uj   3X-  \  2   /+  pQ   3X.J" 


r\j  r>^j     r\j 


buc,  =  v 


3X.3X. 

«J     J 


(1.6) 


where  the  buoyancy  (ft)  is  given  by 


Bf  = 


p0  -  p 


(1.7) 


Assuming  horizontal  homogeneity  of  mean  variables,  where  the  horizontal  mean 

is  defined  by 

X  X^ 

2  r   2 


(?) 


(z,t)  =  lim—  /    /    f (x,y,z,t)  dxdy 

x~° x2  U  h 


2         2 
and  separating  all  variables  into  mean  and  fluctuating  components  gives 


W  - 


Ug  +  u(z,t)  +  u(x,y,z,t) 

Vg  +  V(z,t)  +  v(x,y,z,t) 

0  +  w(x,y,z,t) 


T(z,t)    +  e(x,y,z,t) 


•  Pg 


+  p(x,y,z,t) 


B(z,t)    +  b(x,y,z,t) 


Strictly  speaking,  the  total  fluctuating  part  of  each  variable  includes 
a  component  directly  attributable  to  surface  and  internal  wave  motion.  This 


8 


component  has  been  removed  from  the  fields  depicted  above.  The  virtually 
irrotational  wave  velocities  are  assumed  noninteractive  except  as  an  external 
source  or  sink  of  turbulent  energy  due  to  the  net  contributions  of  breaking 
surface  waves  and  radiating  internal  waves. 

In  practice,  a  time  mean  is  employed  in  data  analysis  rather  than  even 
the  horizontal  mean.  Ideally,  the  averaging  time  should  be  short  compared 
with  the  time  necessary  for  significant  changes  in  the  mean  fields  but  long 
compared  with  the  integral  time  scale  of  the  turbulence. 


At 


+ 

/ 

./At 


At 

2 


B(z,t)  =   (t>)  ~  75T  I  b  dt,  for  example. 

fAt 
2 


In  making  the  Boussinesq  approximation,  the  hydrostatic  component  of 
pressure  has  already  been  eliminated.  The  remaining  mean  pressure  is  assumed 
to  be  geostrophic: 

a=  onfVn.  (1.8a) 


ax    -o   "g 

9P 


g  _ 


=  -  Pn  f   U  .  (1.8b) 


dy  ^    g 


Subtracting  the  geostrophic  equations  from  the  total  momentum  equation 
and  dropping  the  negligible  viscous  terms,  the  equations  for  the  mean  momen- 
tum become  in  complex  notation  (for  the  sake  of  brevity) 


|c.i1fC.|a,  (i.9) 


where  C  =  U  +  i  V  and  c  =  u  +  i  v. 

Using  equations  (1.3),  (1.4),  (1.5),  and  (1.7)  and  again  neglecting 
molecular  fluxes,  a  buoyancy  equation  is  formed: 


t^t 


If  +  "i  iir  -  .rH  •  (1-10) 


Equation  (1.11)  is  the  mean  buoyancy  equation. 


3B 

3t 


3bw 

3Z 


poCp 


(1.11) 


The  use  of  a  buoyancy  equation  reflecting  the  combined  conservations  of 
thermal  energy  and  dissolved  material  is  not  only  more  general  than  a  heat 
equation  alone,  but  it  makes  more  obvious  the  coupling  with  the  mechanical 

energy  budget. 


Using  the  decomposition  into  mean  and  fluctuating  parts  and  taking  the 
mean  of  equation  (1.6)  yields  the  mean  equation  for  the  total  mechanical 
energy.  Where  E  =  u-u.  =  u2  +  v2  +  w2, 


/U2  +  V2  .  E_\  ,  d_ 
3t  I   2     2I3Z 
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3Z 


U  uw  +  V  vw 


bw 


3U 


9ui 
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(1.12) 


The  viscous  diffusion  and  viscous  dissipation  of  the  mean  kinetic  energy  are 
negligible  and  have  been  dropped  in  equation  (1.12). 

The  turbulent  kinetic  energy  equation  is  formed  by  subtracting  the 
scalar  produce  of  (U-j)  and  equation  (1.9)  from  equation  (1.12). 


1_  3£  +  3_ 
2  at   3z 


wT^+l 


—  3U  L  —  3V 

uw  —  +  vw  — 

3Z       3Z 


= 

bw 
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e 

(1.13) 


The  budgets  for  the  individual  components  of  turbulent  energy  are 
formed  in  the  same  manner  from  equation  (1.1)  by  setting  i  =  1,2,3  without 

summation: 


1  3u: 


2  3t 


-  uw 


3U 

3Z 


1  3V: 


2  at 


3Z 


3U_  ,  „   — 
3X 


+  ^3  uv  -  n2  uw  "  3" 


p0  3y  a< 


UV  - 


(1.14a) 


(1.14b) 
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1    1^.^.8.  [wl+.ffil  +£-|5+o2ui7.^.  (1.14c) 


2    ?T  "  uw  "  97  \  T  T  £7  /  T  p7  37  T  "2  uvv  "  3 


on 


Notice  that  the  orientation  of  the  horizontal  axes  is  (x)  positive  to  the 
east  and  (y)  positive  to  the  north^  The  rate  of  viscous  dissipation 
e  =  v  3ui9ui/3xj9Xj  behaves  like  (E/te)  and  the  small  size  of  the  dissipati, 
time  constant  (xe  ^  h//T)  compared  with  the  time  scale  of  the  meteorological 
forcing  causes  the  intensity  of  the  turbulence  to  track  along  in  a  quasi- 
steady  state  with  a  continually  changing  net  rate  of  production.  Hence  the 
time  rate  of  change  of  turbulent  kinetic  energy  is  usually  much  smaller  than 
the  other  terms  of  (1.13)  and  thus  may  be  neglected.  Because  viscous  dissi- 
pation of  energy  occurs  primarily  in  the  range  of  wave  numbers  that  exhibit 
local  isotropy  (the  equilibrium  range),  dissipation  (e)  is  divided  equally 
among  the  component  budgets  (1.14  a-c). 

The  second  term  in  equation  (1.13)  is  the  divergence  of  the  turbulent 
flux  of  kinetic  energy.  Over  the  whole  mixed  layer,  it  probably  accounts  for 
a  net  gain  of  energy.  Wind-wave  interactions  at  the  surface  result  in  some 
net  downward  flux,  primarily  from  breaking  surface  waves.  If  the  Brunt- 
Vaisala  frequency  (N)  of  the  adjacent  underlying  stable  water  column  is  suf- 
ficiently large  so  as  to  be  comparable  with  the  frequency  of  the  integral 
scale  of  the  turbulence,  turbulent  energy  may  be  lost  to  the  generation  of 
internal  waves.  One  of  the  most  significant  aspects  of  this  term  is  that 
locally,  at  the  bottom  of  the  layer  during  occasions  of  entrainment,  a  net 
convergence  of  flux  of  energy  is  necessary  to  maintain  the  downward  buoyancy 
flux  for  a  deepening  mixed  layer. 

The  third  term,  the  rate  of  mechanical  production,  is  perhaps  the 
dominant  source  of  turbulent  kinetic  energy.  It  is  the  rate  of  conversion  of 
mean  to  turbulent  kinetic  energy  by  the  turbulent  flux  of  momentum  down- 
gradient. 

The  last  term  on  the  left,  the  buoyancy  flux,  locally  within  the  mixed 
layer  can  be  either  a  sink  or  a  source.  Usually,  however,  the  mixed  region 
is  slightly  stable  overall,  and  this  term  represents  the  rate  of  increase  of 
potential  energy  by  fluxing  buoyancy  downward.  During  instances  of  large 
buoyancy  flux  up  across  the  surface,  this  term  can  become  an  important  source 
as  in  the  case  of  strong  convective  cooling  in  the  autumn. 

The  summation  of  separate  component  equations  yields  (1.13),  but  one 
term  that  is  wery   important  in  mixed  layer  dynamics  sums  to  zero  and  there- 
fore appears  only  in  the  component  budgets  (1.14  a-c).  This  term  is  the 
correlation  between  pressure  and  rate  of  strain,  p/p0  9ua/9xa.  Since  it  sums 
to  zero  by^  continuity,  it  causes  only  a  redistribution  of  energy  among  if2", 
v7,  and  w 


The  individual  turbulent  energy  budgets  also  have  redistribution  terms 
due  to  rotation  of  the  Earth,  but  these  shall  be  neglected  because  of  the 
usually  short  integral  time  scale  in  comparison  with  one  day.  Perhaps  this 
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effect  does  become  significant  in  some  of  the  very   deep  convective  mixed 
layers  that  are  not  limited  in  growth  by  a  permanent  pycnocline. 

Application  of  the  vertical  homogeneity  assumption  to  the  momentum  and 
buoyancy  equations  (1.9)  and  (1.11)  gives  relationships  for  the  turbulent 
fluxes  in  terms  of  the  boundary  conditions  (specified  externally)  alone: 


57  (z)  =  cw  (0)   (!+£)  +  AC  (|)  jj£ 


(1.15) 


bw  (z)  =  bw  (0)    (1  +|)  +| 


dt      PoCp 
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Qdz 


'owp 


/. 


Qdz 


(1.16) 


The  integral  of  (1.16)  over  the  mixed  layer  gives  the  net  buoyant  damping  for 
the  whole  layer. 
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(1.17) 


Integrating  the  turbulent  kinetic  energy  equations  from  z  =  (-h-6)  to 
z  =  0  gives 
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i^(<v2>h)  = 
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±^(<w*>h)  = 


/0     ^+^lf)dz- 
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The  surface  boundary  conditions  are  prescribed  functions  of  time 


-  uw  (0)  =  — 


0 


(1.20a) 


_      -v  (t) 
=  vw  (0)  =  -a 


(1.20b) 


-  bw  (0)  =  g  [a  sw  (0,t)  -  6  ew  (0,t)] 


(1.21) 


Also  to  be  prescribed  in  deriving  the  system  is  the  radiation  absorption 
function,  Q  (z,t). 

The  boundary  conditions  at  the  bottom  of  the  mixed  layer,  z  =  -h,  will, 
on  the  other  hand,  conform  to  the  developing  situation.  To  derive  these 
conditions,  the  equations  (1.9)  and  (.1.11)  are  integrated  over  the  entrap- 
ment zone  from  z  =  (-h-6)  to  z  =  (-h): 


6/h  Um  0 
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-h 
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1  im 
5/h  -0   /    f  V  dz  =  0 


8UW  ,      i    uv 

t-—  dz  =  -  uw  (-h) 

d  Z 


Therefore 


notation, 


Also, 


JL 

-  uw  (-h)  =  AU  -nr  ,     and  similarly 


JL 

vw  (-h)  =  AV  -tz  ,  or  in  complex 


cw  (-h)  =  AC  ||~  .  (1.22) 


bw  (-h)  =  AB  ^  (1.23) 


where  AC  =  All  +  iAV  and  AB  are  the  respective  jumps  in  the  values  of  the 
mean  variables  across  the  density  interface  separating  the  mixed  layer  from 
the  nonturbulent  region  below.  The  discontinuity  need  not  be  a  perfect  one 
(o  =  0)  for  the  boundary  conditions  (1.22)  and  (1.23)  to  be  valid.  A  suffi- 
cient condition  is  for  the  fluxes  of  momentum  and  buoyancy  from  the  mixed 
layer  into  the  interface  zone,  resulting  in  a  lowering  of  the  vertical  pos- 
ition of  the  zone,  to  be  much  larger  than  that  portion  of  the  fluxes  con- 
tributing to  changes  in  the  momentum  and  buoyancy  profiles  of  the  moving 
interface  zone  itself. 

Integrating  equations  (1.9)  and  (1.11)  from  (z  =  -h  -  6)  to  (z  =  0) 
provides  a  form  of  the  equations  that  includes  the  effects  of  the  entrainment 
stress  and  entrainment  buoyancy  flux  of  (1.22)  and  (1.23). 

h  ^  +  AC  ^  A  =  -  if<C>h  -  c^  (0)  (1.24) 
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and 


h  d<B>    .    AD  dh   . 

ir +  AB  at A 


=  _§£ 


poCp 


J-h-6 


Q  dz  -  bw  (0) 


(1.25) 


where  the  Heaviside  unit  step  function,  A,  is  dependent  upon  (dh/dt) 


rdh^ 


|dt; 


lfor&>0 


0  for  |<0 


(1.26) 


and  <  >  denotes  a  vertical  mean  for  the  mixed  layer: 

0 


.+6  J 


<C>  =  y^  /   C  dz>  etc* 
-h-6 


1.4  Course  of  Action  in  Attacking  the  Problem 

To  lay  a  foundation  and  present  a  perspective  of  the  problem  at  hand, 
the  literature  treating  models  of  the  surface  mixed  layer  is  reviewed  rela- 
tive to  the  basic  principles  and  general  equations  laid  down  in  the  previous 
sections.  This  approach  organizes  the  historical  work  in  terms  of  the  funda- 
mentals, and  it  provides  the  stepping  stones  for  the  development  of  this 
research. 

The  turbulent  kinetic  energy  budget  is  examined  closely.  The  role  of 
the  previously  neglected  redistribution  terms  is  assessed.  All  of  the  terms 
are  modeled  by  use  of  mean-turbulent-field  techniques,  permitting  the  even- 
tual implicit  solution  for  the  turbulent  energy  content  of  the  mixed  layer. 

The  final  preparatory  work  needed  to  complete  the  model  is  treated  in  a 
chapter  on  entrainment.  This  includes  the  derivation  of  an  equation  relating 
the  rate  of  entrainment  (dh/dt)  to  the  other  variables. 

The  final  numerical  method  of  solution  of  the  nonstationary,  non-linear 
set  of  equations  permits  the  solution  of  hypothetical  cases  as  well  as  the 
simulation  of  field  observations.  The  numerical  model  requires  as  input  the 
initial  conditions  of  density  and  current  and  the  surface  boundary  flux  of 
buoyancy  (heat  and/or  salinity)  and  surface  wind  stress  as  functions  of  time. 
Model  outputs  include  the  mean  density  profile,  the  turbulent  kinetic  energy, 
and  the  mixed-layer  depth,  all  as  functions  of  time. 
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2.  REVIEW  OF  THE  LITERATURE 

2.1  Ekman  Depth  of  Frictional  Resistance 

V.  Walfrid  Ekman  (1905)  originated  the  concept  of  a  "depth  of  frictional 
resistance"  for  the  upper  section  of  a  wind-stressed  ocean.  This  depth  (d) 
comes  from  the  mathematical  solution  to  the  steady  state  horizontal  momentum 
equation,  (1.9),  in  which  the  Reynolds  stress  is  related  to  the  mean  shear  by 
a  constant  eddy  viscosity  (K). 
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^  =  -  ifC  - 
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(1.9) 


where 
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—   v    3C 
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If  the  boundary  conditions  are 


-  cw  (0)  =  u|  +  i  0 


and 


-  cw  (-»)  =  0  , 


the  solution  is 
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sin  (  d~   4  } 


/  ttZ  ,  ir  \ 
1  cos  (  d~   4  } 


TTZ 
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(2.2) 


where 


/  2K 


(2.3) 


At  the  depth  (z  =  -  d)  the  direction  of  the  flow  is  opposite  to  the  surface 
current  and  the  magnitude  has  been  reduced  to  (e-7T)  times  the  surface  magni- 
tude. 

Classical  thought  has  suggested  that  any  surface  layer  mixed  by  the 
action  of  the  wind  should  have  a  depth  that  is  of  the  same  order  as  (d). 
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Using  a  quadratic  law  relating  windspeed  (W)  to  surface  stress 

P0uI  =  CDPaW2  (2.4) 

and  an  empirical  relation  from  observations, 

|C(o)|-SJfl2LWf  (25) 

/sin  (j> 
Ekman  derived  a  formula  for  (d)  as  a  function  of  wind  speed  and  latitude  (<j>) 


d  =  7.5  sec"1         W       .  (2.6) 

/sin  <f> 


2.2.  Rossby  Number 

Rossby  and  Montgomery  (1935)  pointed  out  that  the  depth  (h)  of  a  surface 
drift  current  layer  and  Ekman' s  depth  of  frictional  resistance  (d)  are  not 
necessarily  comparable:  the  depth  (h)  has  a  definite  physical  meaning,  but 
(d)  designates  only  the  theoretical  rate  of  exponential  decay  for  a  system 
obeying  (2.1), 

Rossby  and  Montgomery  derived  the  formula  (h  «  W/sin  <f>)  or,  equivalently 

h  -  ujf  (2.7) 

where  the  constant  of  proportionality  is  the  Rossby  number,  R0  =  u*/hf.  They 
then  presented  measurements  demonstrating  the  greater  validity  of  (2.7)  in 
comparison  with  (2.6.) 

It  should  be  recognized,  however,  that  Ekman' s  result  differs  from  that 
of  Rossby  and  Montgomery  only  because  of  the  use  of  the  relation  (2.5.)  If 
instead  of  applying  this  empirical  constraint,  the  eddy  viscosity  (K)  is 
modeled  in  terms  of  likely  turbulent  length  and  velocity  scales  (K  *v*  u*h)  and 
is  assumed  to  be  constant  with  depth,  then  (d  =  2ir2u*/f)  and  the  quadratic 
stress  law  gives 
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So,  in  spite  of  suggesting  a  distinction  between  the  mixed  layer  depth 
and  the  depth  of  fractional  resistance,  Rossby  and  Montgomery  obtained  a 
result  perhaps  more  comparable  with  that  of  Ekman  than  with  the  real  situ- 
ation. This  was  the  case  because  both  derivations  considered  only  the  momen- 
tum budget,  neglecting  the  effects  of  buoyancy  and  mechanical  energy  upon  the 
vertical  fluxes  of  buoyancy  and  momentum  within  this  turbulent  oceanic 
boundary  layer. 


2.3  Eddy  Transfer  Coefficients  in  a  Steady-State  Problem 

Munk  and  Anderson  (1948)  first  combined  the  two  problems  of  density 
structure  and  current  structure  into  a  unified  theory  on  a  steady-state 
thermocline.  Like  Ekman,  they  proposed  an  eddy  viscosity  (Km)  plus  an  eddy 
conductivity  (Kh),  but  these  parameters  were  made  variable  with  the  local 
gradient  Richardson  number  (Ri). 

Ri  =  MGSL  .  (2.9) 

(8U/8Z)2 

This  model  therefore  included  some  of  the  effects  of  the  turbulent  energy 

budget. 

Vh  -  (i  +  c°  „  Ri)  Vh  •  (2-10) 

'    v     m,H 

This  function  for  the  eddy  viscosity  and  eddy  conductivity  was  chosen  because 
of  its  asymptotic  behavior  for  small  and  large  values  of  (Ri): 

1 im    K  u  =  KO,  coefficient  for  no  density  gradient 

Ri  -*  0     ™'H 

1 im    K  n  =  0,  for  extreme  stability. 

Ri  -v  °°     ' 

Because  Munk  and  Anderson  assumed  steady  state  and  did  not  recognize  the 
presence  of  a  sharp  interface  marking  a  boundary  between  the  fully-turbulent 
mixed  layer  and  the  essentially  quiescent  stable  region  below,  their  results 
still  resembled  Ekman' s  original  solution  more  than  they  do  the  physical 
real ity. 


2.4  Obukhov  Length  Scale 

The  more  recent  efforts  in  modeling  the  oceanic  mixed  layer  started  with 
a  one-dimensional  steady-state  study  by  Kitaigorodsky  (1960).  Assuming  that 
the  ocean  surface  mixed  layer  was  analogous  to  the  constant-flux  atmospheric 
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surface  layer  (as  in  Businger  et  al.,  1971),  Kitaigorodsky  concluded  by 
dimensional  analysis  that  the  mixed  layer  depth  must  be  proportional  to  the 
Obukhov  length,  L. 

For  Kitaigorodsky' s  assumptions,  the  momentum  and  buoyancy  equations, 
(1.9)  and  (1.11),  reduce  to 


8Z 

and  8bw 


'  =  0  (2.11) 

=  0  .  (2.12) 


3Z 

The  radiation  absorption,  Q,  was  assumed  to  be  confined  to  the  immediate 
surface  layers.  Taking  the  x-direction  to  be  in  the  direction  of  the  wind, 
the  solutions  to  equations  (2.11)  and  (2.12)  are 

-  cw  =  constant  =  u| 

and  -  bw  =  constant  =  u*b* 

where  (uf)  and  (u*b*)  are  the  downward  surface  fluxes  of  momentum  and  buoy- 
ancy. 

If  the  depth  of  the  surface  mixed  layer  (h)  is  dependent  only  upon  the 
two  parameters  (u*)  and  (b*),  then 

u,2 


(hsg|ew(0)|'    (h  b  ' 


h  =  H*  L  (2.13) 

where  L  =  uj/b*  (2.14) 

and  H*  is  a  constant  of  proportionality. 

If  the  coriolis  force  is  a  significant  component  of  the  mean  momentum 
budget,  then  (2.11)  is  replaced  by 


9cw 

8Z 


=  -  ifC  .  (2.15) 
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Adding  the  coriolis  parameter  (f 
makes  H*  variable. 


=  2Q   sin  <j>)  to  the  dimensional  analysis 


j*  =  — 


(2.16) 


with  a  second  dimensionless  product  (b*/u*f). 

Using  data  from  the  NORPAC  expedition,  Kitaigorodsky  found  that  equation 
(2.13)  with  a  constant  H*  was  insufficient  for  cases  varying  over  more  than 
twenty  degrees  of  latitude  (see  Fig.  3). 

Using  Kitaigorodsky1 s  data,  it  should  be  noticed  that  the  Rossby  number 
(Ro)  based  upon  the  layer  depth  and  (u*), 


D    fh 


(2.17) 


is  less  variant  than  H*  =  h/L  for  the  same  data.     This  can  be  seen   in  Figures 
3  and  4. 


\ 


\ 


Dota 

\  Ro  dependence:  h  a  W/f 

\  h  a  L,  Obukhov  length: 

haW3/(/3gQ) 
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both  Ro~l   and  H*    {data  {^ftom  ICLtaigosi- 
odUky,    1960). 
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There  is  a  basic  flaw  in  this  model  arising  from  the  assumption  that  the 
ocean  mixed  layer  is  analogous  to  the  "constant-flux"  atmospheric  surface 
layer.  The  ocean,  in  all  probability,  does  have  a  surface  layer  over  which 
fluxes  are  approximately  constant  and  a  quasi-steady  state  does  exist.  How- 
ever, the  depth  of  such  a  layer  would  be  limited  to  at  most  the  upper  ten 
percent  of  the  vertical  extent  of  the  entire  mixed  layer.  With  a  constant 
heat  flux  at  the  surface,  the  mixed  layer  temperature  and  depth  cannot  both 
remain  unchanged. 


2.5  Compensation  Depth  for  Shortwave  Radiation 

Kraus  and  Rooth  (1961)  also  conceived  a  steady-state  model  based  pri- 
marily upon  the  buoyancy  equation.  In  their  model  however,  steady  state  was 
made  feasible  by  balancing  the  short  wave  radiation  input  Q(z)  with  a  net 
surface  heat  loss,  p0Cp9w(0)>0,  by  means  of  evaporation,  conduction,  and 
infrared  radiation.  A  compensation  depth  (hc)  is  the  depth  at  which  a  bal- 
ance is  struck  between  the  surface  heat  loss,  p  Cpew(0),  and  the  total  radi- 
ation absorbed  in  the  layer  above, 


0 
Q  d  z 
h 


f, 

c 
r 

If  Q  =  y  Qo  eYZ 

and  if  buoyancy  is  a  function  of  temperature  alone,  then  (1.11)  becomes  in 
steady  state 


8  6W 


8  2   pocp 


+  -L  .  (2.19) 


The  depth  (z  =  -  hc)  is  that  level  at  which  the  turbulent  flux  ew  goes  to 
zero,  with  a  stable  temperature  profile  below  and  an  unstable  one  above. 
Therefore  in  the  region  (0  >  z  >  -  hc)  turbulent  kinetic  energy  can  be  con- 
vectively  produced  since  here 

bw  =  pg  ew  >  0  . 
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Integrating  (2.19)  from  (z  =  -  hc )  to  (z  =  0)  and  using  (2.18)  gives 


hc  =  y"1  In 


1 


1  -  PCpew(0)/Q0 


(2.20) 


Figure  5  is  a  schematic  portrayal  of  the  Kraus  and  Rooth  concept.   If 


ew  =  Kn  —  ,  KH  ^  constant  in  z  , 


(2.21) 


then  (hc)  is  also  the  depth  of  maximum  temperature. 

Using  the  eddy  conductivity  closure  posed  by  (2.21),  Kraus  and  Rooth 
examined  the  structure  of  the  temperature  field  from  the  surface  down  through 
the  mixed  layer  and  across  the  interface,  and  its  variation  with  changes  in 
(Q)  and  the  boundary  conditions.  Their  solutions  are  only  qualitatively 
useful  because  not  only  is  (K^,)  unknown,  but  also  it  is  assumed  to  be  con- 
stant—even across  the  density  interface. 


-hc 


c  =  Y  'An  ( 


h,  = 


Qo^ocpew(0) 


) 


h   >  h. 


/    /    / 
/    /    / 

/    /    /  z 

/    /    / 


solar  radiation 
through  surface 

00>Poc.J^{0) 


conduction,  evap- 
oration and  infra- 
red radiation 

ewio)  >  o 


9w 


6w(z)  =  0w(O)-      °  (l-eyZ) 
rogp 


Vi.QuJm  5.      E{)^cX  of,  the  domp2.n6cutlon  depth,   hc,  cutum-ing  &tmdy  Atate.  and 

Q  =  YQ0eYZ. 
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With  regard  to  the  depth  of  the  mixed  layer  (h),  all  they  could  really 
say  was  it  must  reach  some  depth  greater  than  (hc).  "Dependent  upon  the 
intensity  of  the  convective  and  wind  drive  turbulence,  this  convective  regime 
may  penetrate  more  or  less  deeply  beyond  the  level  (hc)." 

They  visualized  a  steady-state  layer  depth  as  being  possible  only  by 
requiring  upwelling  (W>o)  of  sufficient  magnitude  to  maintain  the  vertical 
position  of  the  entraining  interface. 


2.6  Prototype  Turbulent-Energy-Budget  Model:  Kraus  and  Turner 

Recognizing  the  limitations  in  application  of  the  Kraus  and  Rooth 
model--no  provision  for  a  possible  downward  surface  heat  flux,  no  account  of 
mechanical  production  of  turbulent  kinetic  energy,  and  the  steady-state 
constraint--Kraus  and  Turner  (1967)  further  improved  and  generalized  this 
kind  of  one-dimensional  model.  Their  model  was  the  first  instance  in  which 
it  was  recognized  that  the  budgets  for  thermal  and  mechanical  energies  could 
be  considered  separately.  This  is  valid  because  the  dissipative  rate  of 
heating  (p0/J»e)  is  several  orders  of  magnitude  smaller  than  either  Q(z)  or 
| p 0Cp  3ew/3z|.  Therefore,  their  model  consisted  of  two  separate  equations-- 
the  neat  equation  and  a  mechanical  energy  equation  in  which  the  net  effect  of 
the  work  of  the  wind  on  the  sea  surface  and  the  viscous  dissipation  within 
the  mixed  layer  are  parameterized.  This  use  of  a  mechanical  energy  equation 
together  with  the  buoyancy  (heat)  equation  and  the  boundary  condition  (1.23), 

-  bw~  (-h)  =  AB  ^  ,  (1.23) 

gave  for  the  first  time  a  closed  set  of  equations  whose  solution  provided 
h(t). 

If  buoyancy  is  a  function  of  temperature  alone  and  constant  in  the  mixed 
layer,  equation  (1.25)  is  applicable.  Using  the  radiation  absorption  func- 
tion (2.18)  it  becomes 

h  4#  +  AT  £r  A  =  -^-  (1  -  e"Yh)  -  ew  (0)  .         (2.22) 
dt      dt     p0C 

The  turbulent  kinetic  energy  equation,  integrated  from  the  top  of  the  en- 
trapment zone  (z  =  -  h)  to  the  surface  is 

f°    — 
G  -  D  =  -3g   /    ew  dz  (2.23) 
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where 


r°- 


w 


3U. 


dz 


is  the  total  rate  of  mechanical  production  and 

0 


J-h 


D  =   /    £dz 


is  the  total  rate  of  dissipation  within  the  mixed  layer.  Neither  of  these 
parameters  in  the  Kraus-Turner  model  is  an  implicit  variable,  and  each  must 
be  specified  externally.  The  heat  equation  that  leads  to  (2.22)  is  also  used 
to  eliminate  ew  from  (2.23)  by  integrating  between  z  and  0: 


m 


96W 
3Z* 


P0Cp 


) 


dz'  =  0 


or 


ew  (z)  =p-f  (1  -  eYZ)  +  ew  (o)  -z^f 


(2.24) 


Integrating  ew  (z)  as  prescribed  by  (2.24)  from  (z  =  -  h)  to  (z  =  0)  gives 


/.' 


7—  j    1  u?  d<T> 
ew  dz  =  75-  hz   ■■ 
c  dt 


OP 


-  ew(0) 


h  + 


Yn  C 

rpo  p 


1  -  e 


-Yh 


(2.25) 


Neglecting  e"Yh,  equations  (2.22),  (2.23)  and  (2.25)  can  be  used  to  give 

another  equation  (2.26),  which  together  with  j_2.22)  constitutes  a  closed 

system  in  h(t)  and  <T>  (t)  where  G,  D,  Q0  and  ew(0)  are  prescribed  functions 
of  time  at  most. 


h2  d<T>   yT.  dh  .    G  -  D 
T  -3T  +   ATh  dt  A  =  -W  + 


PoV 


(2.26) 


An  important  contribution  by  Kraus  and  Turner  was  the  conceptualization 
of  a  model  for  which  a  stationary  or  even  retreating  mixed  layer  depth  is 
possible.  In  such  a  case  where 
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#*»■  A^  =  0 


equations  (2.22)  and  (2.26)  are  still  applicable  because  of  the  presence  of 
the  Heaviside  unit  step  function.  Setting  (h  =  hr)  in  the  case  of  a  retreat- 
ing or  steady  mixed  layer  depth,  these  equations  reduce  to  (2.22a)  and 
(2.26a). 


h2    - 

nr  d<T>   G  -  D 


+  — r-  .  (2.26a) 


dt    eg    p0cpY 


Neglecting  short-wave  radiation  that  escapes  the  mixed  region  and  eliminating 
(d<T>/dt)  between  (2.22a)  and  (2.26a)  gives 

h  =  2(G  "_^ .  (2.27) 

eg  (Q0  -  ew(0)) 

Whenever  the  surface  boundary  conditions  and/or  solar  radiation  adjust 
to  make  (hr  <  h)  the  mixed  layer  will  "retreat."  Of  course,  the  region  does 
not  unmix,  in  accordance  with  the  second  law  of  thermodynamics.  The  net  rate 
of  production  of  turbulent  kinetic  energy,  G-D,  is  insufficient  to  balance 
the  rate  of  increase  of  potential  energy, 

r°       - 

I      -  eg  ew  dz, 
J-h 


required  to  mix  the  region  all  the  way  to  the_  density  interface.  Conse- 
quently, as  the  region  warms  (notice  that  d<T>/dt  can  be  only  positive  at 
this  time),  a  new  density  interface  is  established  at  z  =  -  hr. 

Kraus  and  Turner  model  G  in  terms  of  the  friction  velocity: 

G  =  u*3  .  (2.28) 

Not  knowing  the  importance  of  the  viscous  rate  of  dissipation, 
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D  =  /   £dz, 
J-h 

they  simply  neglect  it.  They  find,  however,  that  their  model  predicts  a  too- 
large  value  for  (h),  and  that  the  dissipation  should  possibly  be  included. 

Turner  (1969)  felt  that  there  was  a  problem  with  (2.28)  as  well.  From 
examination  of  observations  of  sudden  wind  speed  increase  and  ensuing  mixed 
layer  deepening,  he  deduced  that  "a  substantial  fraction  of  the  part  of  the 
work  done  by  the  wind  which  goes  into  the  drift  current  is  eventually  used 
to  deepen  the  surface  layer."  This  statement  reflects  the  need  for  a  more 
comprehensive  model,  particularly  for  unsteady  situations.  Such  a  model 
should  reflect  the  input  of  energy  into  the  mean  velocity  profile  and  the 
time  delay  needed  to  shift  some  of  this  energy  to  turbulence. 

Again  with  regard  to  the  Kraus  and  Turner  model,  the  setting  of  (D  =  0) 

so  that 


0 

-h 


u 
-eg  f       ew  dz  =  G  (2.29) 


places  an  unrealistic  constraint  upon  the  buoyancy  term:  it  becomes  depen- 
dent only  upon  the  mechanical  production.  The  error  in  this  is  most  obvious 
when  there  is  strong  surface  cooling  and 

(-  /   ew  dz) 
J-h 

is  less  than  zero. 

Kraus  and  Turner  also  neglect  the  effect  of  entrainment  in  their  turbu- 
lent kinetic  energy  budget,  equation  (2.23). 

In  spite  of  these  deficiencies,  this  model  was  a  big  step  in  the  right 
direction  in  its  consideration  of  the  turbulent  energy  budget  in  recognition 
of  the  energy  source  for  mixing  and  entraining. 


2.7  Adding  Dissipation 

MiropoTskiy  (1970)  and  later  Denman  (1973a)  assume  that  dissipation  is 
a  fixed  fraction  of  the  shear  production,  equation  (2.30). 

D  =  6'  G  (2.30) 
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where  6'  is  an  empirical  parameter  to  be  determined  from  observations 
gives,  instead  of  (2.29),  the  equation  (2.31): 


This 


-  eg 


/ 


0 


ew  dz  =  (1  -  6')  G 


(2.31) 


This  does  not  solve  the  dissipation  problem  because  6"  cannot  have  a  constant 
value.  In  essence,  dissipation  must  be  allowed  to  adjust  to  the  total  situ- 
ation as  it  evolves. 

Miropol'skiy  also  assumed  that  G°=u|,  but  a  variation  of  an  exercise  he 
uses  to  deduce  this  demonstrates  perhaps  the  major  source  of  error  in  a  model 
like  (2.28). 


G  =  - 


L  (uwfl 


+  vw  —   dz 


If  3uw/3z  =  -  fV  and  3vw/9z  =  fU,  then 


In  general ,  however. 


-  uw  V  -  uw  V 


-h 


3uw  _   ...   all     , 
fV  -  w  ,  and 


9VW    .c,,    3  V 

JT=   fU  "  at  '  giving 


^  U: 


G  = 


uw  U  -  vw  V 


0 
-h 


(2.32) 


and  thus  indicating  the  importance  of  the  mean  kinetic  energy  and  its  distri- 
bution within  the  mixed  layer.  As  will  be  shown,  even  if  the  wind  is  steady 
for  long  periods  of  time,  the  mean  kinetic  energy  can  change  markedly  on  a 
time  scale  corresponding  to  the  inertia!  period. 

Most  recently,  the  trend  in  the  literature  has  been  to  model  the  mixed 
layer  as  a  vertically  homogeneous  moving  slab  with  density  and  velocity 
discontinuities  at  the  entraining  interface.  Geisler  and  Kraus  (1969)  were 
the  first  to  use  the  slab  approach  in  their  model  of  the  atmospheric  boundary 
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layer.  This  problem  is  almost  completely  analogous  to  the  oceanic  mixed 
layer  problem  except  that  the  atmospheric  boundary  layer  is  driven  by  a 
horizontal  pressure  gradient  rather  than  by  a  surface  stress.  A  rigid  sur- 
face boundary  rather  than  a  free  surface  also  results  in  a  subtle  but  impor- 
tant difference.  Nevertheless,  the  basic  setups  of  the  two  problems  are  the 
same. 

Simultaneous  calculation  of  the  mean  velocity  together  with  (h)  and  (T) 
permits  the  implicit  calculation  of  the  mechanical  production  of  turbulent 
energy,  an  improvement  upon  the  previous  (u|)  method. 

The  equations  (2.33)  and  (2.34),  reflecting  the  conservations  of  mean 
momentum  and  mean  buoyancy  are  essentially  the  same  as  (1.24)  and  (1.25). 
The  only  real  difference  is  that  Geisler  and  Kraus  assume  a  prescribed  mean 
subsidence  (analogous  to  ocean  upwelling)  in  their  atmospheric  boundary 
layer.  This  non-zero  vertical  velocity  (W)  can  result  in  a  stationary  mixed 
layer  depth  even  when  entrainment  is  occurring.  This  then  is  a  different 
mechanism  than  that  developed  by  Kraus  and  Rooth  (1961)  for  obtaining  a 
constant  (h). 

h  ^^  +  AC  (~  -   W)  =  -  ifh  (<C>  -  Cg)  -  cw  (0)  .         (2.33) 


h  ^T+  AT  (dT-  W)  =  "  ™  (0)   •  (2'34) 

In  (2.33),  (-  ifCg)  is  the  kinematic  geostroplrk  pressure  gradient,  the 
source  of  momentum.  The  kinematic  surface  stress,  cw  (0)  is  a  momentum  sink 
in  this  case. 

Geisler  and  Kraus  seem  to  avoid  the  problem  of  dealing  with  the  viscous 
dissipation  of  turbulent  kinetic  energy  by  prescribing  a  fixed  value  for  the 
integrated  flux  Richardson  number,  Rfj. 

.h  _ 
bw  dz 


/„ 


RfT  =  — s  ~i  —x—  ■=  n,  constant  .        (2.35) 


f  h  /  UW  3U  "  UW  9V\ 
JQ       \  3Z       ZZ) 


dz 


This,  however,  is  equivalent  to  Miropol 'skiy's  method.  From  (2.31)  and 
(2.35), 

n  =  1  -  5'  .  (2.36) 

Since  (n)  is  always  a  fixed  positive  number,  and  with  no  radiation 
absorption  in  the  model,  the  buoyancy  flux  (heat  flux)  can  only  be  downward. 
Therefore,  this  model  like  that  of  Ki taigorodsky  is  restricted  in  application 
to  only  those  cases  where  the  mixed  layer  is  stable  throughout. 
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2.8  Role  of  Mean  Kinetic  Energy 

Pollard,  Rhines  and  Thompson  (1973)  apply  the  slab  approach  to  the 
oceanic  mixed  layer,  but  they  complete  the  entrainment  problem  with  a  differ- 
ent mechanical  energy  requirement.  The  time  rate  of  change  of  the  total  mean 
mechanical  energy,  potential  plus  kinetic,  of  the  slab  is  set  equal  to  the 
rate  of  work  by  the  wind  on  the  mean  flow: 
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(2.37) 


This  is  the  same  as  (1.12)  integrated  from  (z  =  -  h  -  6)  to  the  surface  if 
viscous  dissipation  is  neglected  and  (U)  and  (V)  are  constant  in  (z)  within 
the  mixed  layer.  The  turbulent  buoyancy  flux  (bw)  gives  the  potential 
energy  change  using  equation  (1.11): 
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Neglecting  the  time  rate  of  change  of  the  turbulent  kinetic  energy, 
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If  radiative  heating  is  ignored,  (2.37)  reduces  to 


* 
Ri  = 


hAB 


=  1 


U2  +  V2 
Pollard,  Rhines  and  Thompson  assume  that  (2.37)  applies  as  long  as 


(2.37a) 
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is  positive.  As  soon  as  this  rate  of  work  by  the  wind  becomes  negative, 
"energy  flow  to  increase  (h)  ceases  and  since  the  water  cannot  unmix,  (h) 
must  be  constant...".  Therefore,  mixed  layer  deepening  would  occur  only  up 
until  one-half  of  an  inertial  period  following  the_  onset  of  a  steady  wind 


stress.  This  result  is  demonstrated  by  setting  cw  (0)  = 
AC  =  <C>  in  equation  (1.24),  giving 


ul  and 


*hf>hi+  if  (<c>h) 


=   U: 
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The  particular  solution  is  <U>  +  i  <V>  =  (iu^/fh)(e-1 ' t   -1).  Hence  x-fr  (z 
=  0)  goes  to  zero  when  <U>  vanishes  at  time  t  =  -n/f. 

In  this  model,  energy  for  entrainment  is  derived  directly  from  the  mean 
flow.  A  separate  budget  for  turbulent  energy  is  not  considered.  That  is  to 
say,  the  intensity  of  the  turbulence  is  not  considered  to  have  an  active  role 
in  the  mechanism  of  entrainment  in  this  model. 

In  a  three-layer  model  of  the  ocean  mixed  layer,  Niiler  (1974)  combined 
elements  of  Pollard  et  al .  and  Kraus  and  Turner  (1967)  in  that  both  turbulent 
kinetic  energy  and  the  mean  kinetic  energy  are  considered  to  be  important  in 
the  mechanism  that  determines  rate  of  entrainment.  Figure  6  is  a  diagram  of 
the  vertical  temperature  and  velocity  structures  in  this  model.  The  turbu- 
lently  active  region  was  divided  into  three  subregions:  (i)  a  constant-flux 
surface  layer,  (ii)  the  major  part  of  the  whole  region,  and  (iii)  the  en- 
trainment zone,  lying  just  above  a  "quiescent  abyss." 

The  momentum  and  buoyancy  equations  used  by  Niiler  are  virtually  the 
same  as  the  "slab"  equations  (1.24)  and  (1.25)  if  B  =  B(T)  only. 
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h  d<C>   ,    .    .   r  dh       .£U     r        —  /n\       F 

—jT—  +   A   A   C  -rr  =    lfh    <C>  -   CW    (0) . 

at  at  p 


(2.38) 


h  d<T>  _   .    .   T  dh  _      —  /nx 

ir_AATdt      9w(0) 


(2.39) 


Radiation  absorption  was  assumed  to  occur  entirely  in  the  uppermost 
layer  and  was  therefore  included  in  ew  (0).  The  additional  term  (F/p0)  is 
a  "damping"  force  for  inertial  motions 
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(2.40) 


within  the  mixed  layer  and  presumably  is  related  to  the  "'dissipation'  of 
mean  motions  as  well  as  the  radiation  flux  of  momentum  from  the  bottom  of  the 
mixed  layer."  Pollard  and  Millard  (1970)  considered  such  a  term  as  well,  but 
one  that  was  linear  and  thus  resulting  in  an  exponential  damping.  The  time 
constant  (tj)  ranged  from  four  to  twenty-five  inertial  periods,  depending 
upon  the  size  of  the  inertial  circle  relative  to  the  horizontal  scale  of  the 
forcing  wind  system. 
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If  the  mean  velocity  below  the  density  interface  is  zero  and  the  mean 
temperature  below  the  interface  is  given  by  T(z  <  -  h)  =  rTz,  then  the  fluxes 
at  (z  =  -  h)  reduce  to 
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.r   dh    r     dh 
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(2.41) 


and 


e^  (-h)  =  a  T  ^  =  (<T>  +  rTh)  {£  . 


(2.42) 


Since  T  =  <T>  and  C  =  <C>  for  the  bulk  of  the  mixed  layer,  the  turbulent 
fluxes  are  linear  functions  in  (z): 
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(2.43) 
<T>  +  rTh)  ^J.   (2.44) 


The  relative  importance  of  the  terms  of  the  turbulent  kinetic  energy  budget 
(1.13),  was  hypothesized  to  vary  with  subregion. 
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The  rate  of  mechanical  production,  -un-w  aU-j/az,  was  assumed  to  be  non-zero 
only  in  the  surface  and  entrainment  subregions.  In  these  two  regions  of 
small  vertical  extent,  the  buoyant  damping  was  considered  by  Niiler  to  be 
insiqnificant.  Thus  the  balance  is  between  mechanical  production,  turbulent 
and  pressure  diffusion,  and  dissipation,  equation  (2.45a): 
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(2.45a) 


In  this  model  aU/az  =  aV/az  =  0  within  the  central  part  of  the  mixed  region, 
and  hence  mechanical  production  is  necessarily  assumed  to  be  zero.  There- 
fore, here  the  buoyant  damping  and  viscous  dissipation  are  balanced  by 
diffusion  from  the  two  adjacent  production  layers,  giving  equation  (2.45b). 
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Integrating  equations  (2.45a)  and  (2.45b)  vertically  and  combining  them  to 
eliminate  w(E/2  +  p/pq)  at  (z  =  -  d')  and  (z  =  -  h)  gives 


^1-  /d°(-f+-§) 


J 


aV  , 
vw-\dz  • 


dz 


+    1<C>I2     dh  = 


dt 


/' 

^-h 


bw  dz   . 


(2.46) 


In  the  manner  of  Kraus  and  Turner,  Niiler  parameterized  the  sum  of  the 
first  three  terms  of  (2.46)  in  terms  of  the  surface  stress. 
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Integrating   (2.42)  to  give 
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bw  dz  =  -  3g 
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and  using  (2.47),  equation  (2.46)  becomes 
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The  system  of  equations  (2.38),  (2.39),  and  (2.48)  is  a  closed  set  of  equa- 
tions in  the  three  unknowns:  -  h,  <T>  and  <C>. 

This  model  more  resembles  that  of  Kraus  and  Turner  than  it  does  that 
of  Pollard  et  al  because  of  the  utilization  of  a  parameterized  turbulent 
kinetic  energy  equation  rather  than  a  total  mechanical  energy  equation. 
The  primary  difference  is  the  presence  of  the  entrainment  production  term, 
|<C>|2/2«dh/dt,  in  (2.48)  which  necessitates  the  additional  equation,  the 
integrated  momentum  equation  (2.38). 

One  aspect  of  important  consequence  in  Niiler's  model  manifests  the  need 
for  an  even  more  comprehensive  term-by-term  modeling  of  the  turbulent  kinetic 
energy  equation.  This  is  the  fact  that  the  turbulent  kinetic  energy  produced 
by  entrainment,  |<C>2/2«dh/dt,  must  go  entirely  toward  increasing  the  poten- 
tial energy.   Because  of  the  parameterization,  (2.47),  dissipation  is  not 
permitted  to  adjust  to  include  either  the  direct  effect  of  this  particular 
source,  or  the  less  obvious  effects  related  to  the  entrainment  or  lack  of 
it. 

The  only  solution  to  this  predicament  would  seem  to  be  to  model  dissi- 
pation separately  from  any  of  the  source  terms,  allowing  it  to  adjust  to 
total  turbulent  intensity. 


3.  CLOSING  THE  PROBLEM 

The  equations  (1.19a-c),  (1.24),  and  (1.25)  do  not  by  themselves  consti- 
tute a  closed  system  of  equations  for  the  mixed  layer.  Vertical  integration 
over  the  mixed  layer  simplified  the  equations  but  introduced  yet  another 
unknown,  the  mixed  layer  depth  (h). 
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3.1  Net  Viscous  Dissipation  in  the  Mixed  Layer 
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A  dissipation  time  scale  (xe)  is  defined  by  e  =  <E>/xe. 

For  fully-turbulent  geophysical  flows  having  large  Reynolds  numbers, 
viscous  dissipation  of  the  turbulence  occurs  primarily  in  the  small  eddies 
that  are  locally  isotropic.  As  explained  by  Tennekes  and  Lumley  (1972),  an 
inviscid  estimate  of  dissipation  may  be  made  by  taking  the  rate  at  which 
large  eddies  supply  energy  to  small  eddies  (equal  to  the  rate  of  dissipation) 
to  be  proportional  to  the  reciprocal  of  the  time  scale  of  the  large  eddies. 
If  the  time  scale  of  these  large  eddies  is  proportional  to  the  mixed  layer 
depth  divided  by  the  rms  turbulent  velocity  /  <E>,  then  an  integral  model  for 
dissipation  in  the  mixed  layer,  independent  of  viscosity  and  the  small  scales 
is 

•0 
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-h-6 


d  z  =  mi<E>3/2  (3.1) 


wh^ere  (m^  is  a  constant  of  proportionality.  For  those  situations  where 
<E>  =  u|,  equation  (3.1)  is  the  same  as  that  used  by  Miropol'skiy  (1970)  and 
Denman  (1973). 

An  important  concept  in  modeling  dissipation  is  that  of  local  isotropy. 
Turbulent  kinetic  energy  generated  at  the  largest  scale  (^h)  is  transferred 
without  much  additional  production  or  dissipative  loss  through  the  inertial 
subrange  to  the  larger  and  larger  wave  numbers  (smaller  eddies)  by  vortex 
stretching.  Dissipation  is  significant  only  at  the  lower  end  of  this  iner- 
tial subrange. 

Because  there  is  no  preservation  of  the  original  orientation,  the  small 
eddies  of  the  inertial  subrange  are  "locally  isotropic."  Therefore,  dissi- 
pation draws  approximately  equally  from  all  three  turbulent  energy  compon- 
ents. Of  course,  the  existence  of  an  inertial  subrange  is  dependent  upon  a 
large  Reynolds  number,  and  this  is  certainly  the  case  for  the  oceanic  mixed 
layer. 


3.2  Net  Effect  of  Redistribution  of  Turbulent  Energy 

R   =     /     B-  r-2-  dZ  . 
a      /     Pn  3X 

4-6   °   a 
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As  previously  discussed,  Ri  +  R2  +  R3  =  0,  but  Ra  may  be  an  important  source 
of  sink  term  for  the  individual  turbulent  kinetic  energy  budgets. 

Following  the  early  lead  of  Rotta  (1951),  but  in  agreement  with  the 
dominant  term  of  the  rational  closure  technique  of  Lumley  and  Khajeh-Nouri 
(1974), 


R  =  -m2  /e>  (3  <u2- 


:E>) 


(3.2) 


In  addition  to  dimensional  consistency,  the  concept  leading  to  (3.2)  is  that 
of  a  "return  to  isotropy."  In  other  words,  the  correlation  of  pressure  and 
turbulent  rate  of  strain  tends  to  redistribute  energy  equally  among  the  three 
components. 


3.3  Shear  Production 
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where  |6C0|  is  the  "excess"  surface  mean  velocity  in  the  direction  of  the 
wind  stress.  Notice  that  in  this  instance  the  inhomogeneity  of  the  mean 
velocity  field  cannot  be  neglected. 
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where  uj  =  |cw  (0) | . 

If  I  <S C 0 1  1S  proportional  to  (u*),  then  (3.3)  may  be  combined  with  the  param- 
eterized net  input  from  breaking  waves  less  loss  to  radiating  internal 
waves: 
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This  is  basically  in  concurrence  with  the  work  of  Kraus  and  Turner,  Denman, 
and  Niiler. 
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3.4  Need  For  an  Entrainment  Equation 

Using  the  suggested  relationship  (1.24),  (1.25),  and  (3.2),  equations 
(1.14  a-c)   become  (3.5  a,b).  The  two  horizontal  component  equations  have 
been  added  together  to  give  an  equation  for  horizontal  turbulent  kinetic 
energy  <g^>  =  <u->  +  <v*>,  together  with  the  equation  for  vertical  turbulent 
energy  <w->: 
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If  h  (t)  is  unknown,  equations  (3.5a,b)  together  with  (1.24)  and  (1.25) 
are  an  incomplete  system  of  four  equations  in  five  unknowns:  h,  <C>,  <B>, 
<q2> .     An  entrainment  hypothesis  will  provide  the  fifth  equation  needed  to 
close  the  system. 


4.   ENTRAINMENT  HYPOTHESIS 

4.1  Entrainment  in  Earlier  Mixed  Layer  Models 

In  this  study  the  entrainment  velocity,  ue  =  dh/dt,  will  be  modeled 
explicitly  in  terms  of  the  other  free  parameters  of  the  system.  However,  in 
much  of  the  literature  treating  models  for  the  ocean  mixed  layer,  (ue)  is  a 
consequence  of  various  assumed  constraints  placed  upon  the  mechanical  energy 
budget  for  the  layer  as  a  whole. 

In  the  first  really  tenable  model  of  the  mixed  layer  that  was  capable  of 
simulating  a  growing  mixed-layer,  Kraus  and  Turner  (1967)  assumed  that  all 
of  the  turbulent  kinetic  energy  produced  in  the  mixed  layer  goes  to  increase 
the  potential  energy  of  the  system.  After  taking  into  account  any  surface 
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buoyancy  flux,  the  balance  of  the  potential  energy  change  went  to  entrap- 
ment— mixing  the  requisite  amount  of  underlying  denser  water  uniformly 
throughout  the  homogeneous  mixed  region. 

Their  resultant  entrainment  velocity  is 


2  (G  -  D)  +  wb(0)h  +  I-H-)  0/h 

u  = ^P/  (4  1) 

ue  hAB  {*mi) 


where  (G-D)  is  the  net  rate  of  mechanical  production  of  turbulent  kinetic 
energy  minus  the  rate  of  viscous  dissipation  for  the  whole  layer.  The 
surface  buoyancy  flux  wb(0)  may  either  increase  or  decrease  the  entrainment 
velocity,  depending  upon  its  sign.  Again,  not  knowing  how  to  deal  with  the 
dissipation  (D),  Kraus  and  Turner  ignored  it.  They  set  G  =  u|,  so  (4.1) 
becomes 

u  «£ [4. la) 

e      hAB 


where  the  solar  heating  function  (0/)  and  the  surface  buoyancy  flux  have  been 
combined  in  defining  a  buoyancy  flux  scale: 


§221-  wb  (0) 


pnCn 
b*  =  -^ (4.2) 

Geisler  and  Kraus  (1969),  Miropol'skiy  (1970),  Denman  (1973)  and  Niiler 
(1974)  all  assumed  that  a  fixed  fraction  of  the  mechanical  production  of  the 
turbulent  kinetic  energy  would  be  dissipated.  Therefore,  their  resultant 
entrainment  rates  are  the  same  as  (4.1a)  except  that  some  constant  smaller 
than  2.0  would  precede  u|. 

Pollard,  Rhines  and  Thompson  (1973)  did  not  explicitly  consider  the 
turbulent  part  of  the  mechanical  energy  budget,  and  therefore  a  conceptually 
different  relationship  results  from  their  use  of  the  total  mechanical  energy 
equation: 

(hAB  -  |aC|2)  ^jj-  =  u*b*h  .  (4.3) 
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If  (b*  =  0)  this  reduces  to  Ri*  =  1,  or 


h>  AC 


AB 

The  (>)  sign  was  added  to  prevent  the  mixed  layer  from  retreating  when  mean 
kinetic  energy  is  removed  in  the  second  half  to  the  inertial  cycle. 

This  approach  at  first  seems  to  give  a  plausible  result  based  upon  the 
stability  of  the  mean  flow.  The  problem  is  that  the  mixed  layer  is  already 
turbulent,  and  the  system  has  an  excess  of  turbulent  kinetic  energy,  some  of 
which  may  be  available  for  mixing  at  the  interface,  regardless  of  the  value 
of  the  overall  Richardson  number,  Ri*.   It  is  granted  that  Ri*  may  be  con- 
strained to  having  a  value  greater  than  some  critical  value  by  reason  of  a 
mean  flow  instability,  but  the  fact  is  that  in  the  laboratory  and  in  geophys- 
ical cases,  measurements  indicate  that  entrainment  occurs  even  though  Ri*  is 
much  larger  than  one. 

The  entrainment  experiments  with  mean  shear  of  Kato  and  Phillips  (1969) 
and  Moore  and  Long  (1971)  (Figs.  7  and  8)  show  no  direct  relationship  between 
Ri*  and  a  maximum  layer  depth.  More  importantly,  Turner  (1968)  and  others 
have  conducted  experiments  with  growing  mixed  layers  that  had  only  turbulence 

and  no  mean  shear  (Ri*  =  °°). 

All  of  these  laboratory  results  when  examined  together  strongly  imply 
that  the  bulk  Richardson  number  Ri*  is  not  the  parameter  most  relevant  to 
entrainment  rate.  Instead,  a  similar  nondimensional  number,  using  the  tur- 
bulent kinetic  energy  (E)  rather  than  |aC|2,  is  suggested. 

As  shown  by  Niiler  (1974),  the  mean  kinetic  energy  can  influence  entrain- 
ment rate  by  increasing  (E)  in  a  growing  mixed  layer.  This  in  turn  increases 
the  average  position  of  the  bottom  of  an  oceanic  mixed  layer  undergoing  a 
sequence  of  entrainment  and  retreat  due  to  varying  wind  stress  (uj)  and  heat- 
ing/cooling (u*b*)  cycles.  Figure  9  therefore  indicates  some  statistical  de- 
pendence upon  the  vertically-averaged  Richardson  number  and  hence  Ri*  as  well 


4.2  Suggested  Turbulent  Mechanism 

Benjamin  (1963)  shows  that  three  basic  types  of  instabilities  are  pos- 
sible for  a  system  where  a  flexible  solid  is  coupled  with  a  flowing  fluid. 
At  the  interface  between  the  mixed  layer  and  the  denser  water  beneath,  a  so- 
called  class  "A"  instability  will  arise  if 


where  (k)  is  the  wave  number  of  the  interfacial  disturbance, 
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and  Long    [1971)   e.xpeAimznt. 


f|  =  no  e 


ik(x-cxt) 


and  (auO  is  the  total  velocity  change  across  the  interface. 

From  Lamb  (1932),   §232,  a  class  "C"  instability,  the  Kelvin-Helmholtz 
instability,  requires  a  higher  (au) . 


au  >y 


2AB 


However,  the  class  "A"  instability  is  dependent  upon  energy  dissipation  in 
the  lower  fluid,  and  this  is  likely  to  be  small  compared  with  inferred 
rates  of  convergence  of  energy  flux,  -  3/3z  [w(p/p0  +  E/2)]_n,  at  the  inter- 
face. For  geophysical  flows  of  this  type  having  large  Reynolds  and  Pe*clet 
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DEPTH  flvERACCC  RIOARCSCN  NLTBER.    15.8-26-2  METERS 


DEPTH  flVERRCED  RICHARDSON  NUMBER.   26.2-46-2  METERS 
111 


FlguAz   9.  Vzpth-avz/uxgzd  gradient  RichaAdi>on  numbeA,   AzAB/|AC|2, 
v£/l6o4  time,  at  two  dupth  nxmgz&    [buoy  mcaAuAments  comttuy  o& 
Vavld  HatpoAn) .     \Jatuzi>  gucatcn.  than  100  weA&  defined  ai>  equal 
to   TOO.     Vaik&d  kofhizontat  line  at  1.0  izptieAzntA  the.  condition 
^on.  maAgtnat  dynamic  i>tabiJUXy. 


numbers,  the  class  "C"  instability  is  therefore  most  likely  to  be  the  domin- 
ant mechanism  leading  to  observed  rates  of  entrainment. 

The  specific  mechanism  that  is  envisioned  in  the  destabilization  of  the 
interface  and  the  resulting  entrainment  is  a  "local"  Kelvin-Helmholtz  (K-H) 
instability.  The  shear  needed  to  trigger  such  an  instability  is  provided  by 
the  local  turbulent  eddies.  The  mean  shear  contributes  to  the  instability 
but  cannot  in  itself  generate  a  critical  Richardson  number.  The  mean  gradi- 
ent Richardson  number 


Ri 


:    _   9B/3Z 


(3U/3Z)2 


(4.4) 


would  not  be  likely  to  achieve  a  critical  value  because  the  total  instan- 
taneous Richardson  number, 
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R.  =  j^z_  (4>5) 

1    (3U/3Z)2 

is  the  relevant  parameter.  The  minimum  value  of  the  envelope  of  Rij(t)  at 
the  interface  would  determine  the  advent  of  any  appreciable  interface  insta- 
bilities. 

The  onset  of  the  K-H  instability  and  its  exponential  growth  rate  is  pre- 
dicted by  linear  two-dimensional  wave  theory.  As  individual  wave  packets 
achieve  a  significant  amplitude,  the  nonlinear  and  three-dimensional  effects 
of  the  turbulence  field  prevail  by  distorting  the  wave  shapes  and  advecting 
parts  of  the  exposed  cusps  of  denser  water  up  into  the  mixed  layer.  There- 
fore, only  the  initial  stages  of  the  instability  are  strictly  of  the  K-H 
type,  where  the  induced  suction  at  the  crests  of  a  perturbation  wave  on  the 
interface  is  large  enough  to  overcome  the  restoring  buoyancy  force. 


4.3  Relevant  Parameters  and  a  Dimension  Analysis 

Since  the  mechanical  energy  needed  to  continue  this  mixing  process  and 
thus  provide  for  a  significant  entrainment  velocity  must  come  from  the 
turbulent  eddies,  the  rate  of  supply  of  turbulent  kinetic  energy,  -  d/dz 
[w(p/p0  +  E/2)]_.,  just  above  the  interface  should  determine  the  value  of 
(ue)  for  a  given  buoyancy  jump  (aB)  and  mean  velocity  drop  | AC |  across  the 
interface. 

As  Long  (1974)  summarizes  experimental  studies,  the  nondimensional 
entrainment  velocity  W*  is  found  to  depend  upon  the  first  power  of  the 
dimensionless  parameter  E*  based  upon  the  buoyancy  jump  across  the  interface, 
the  depth  of  the  homogeneous  layer  (h),  and  the  intensity  of  the  turbulence 
<E>.  Where  Long  assumes  <E>  a  uf, 


and  W  =  —  aE  .  (4.7) 

The  relationship  (4.7)  together  with  an  integrated  turbulent  kinetic 
energy  equation  to  provide  <E>(t)  plus  the  mean  buoyancy  equation  could  be 
used  to  close  the  problem. 

The  equation  (4.7)  is  appealing  because  it  depends  strongly  upon  that 
which  is  accomplishing  the  erosion  of  the  interface,  the  turbulent  eddies 
having  a  length  scale  (Mi)  and  velocity  scale  (^  /  <E>).  However,  its  use 
based  only  upon  laboratory  evidence  raises  some  questions  about  what  is 
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really  happening  in  the  turbulent  kinetic  energy  budget  and  why  the  mean 
velocity  jump  |aC|  does  not  appear  to  be  important  at  a  first  glance. 

Resorting  to  a  dimensional  analysis  in  which  the  relevant  parameters  are 

h,  u  ,  AB,  |aC|,  and  <E>  gives 

v  (W*,  E*,  Ri*)  =  0  (4.8) 


where 


and 


w*  = 

e 

E*  = 

<E> 
hAB 

Ri* 

hAB 

AC 


Here,  the  physical  significance  of  (h)  is  its  assumed  proportionality  to  the 
length  scale  of  the  turbulent  motion  responsible  for  the  interface  instabil- 
ities. 

Re-writing  equation  (4.8)  so  as  to  solve  for  (ue), 

u  -vG&f-te  ,  -^-)    .  (4.8a) 

e  VhAB        |AC|2/ 

The  experimental  relationship  (4.7)  suggests  that  Ri*  may  be  of  negli- 
gible importance  compared  with  E*.  Other  investigators  besides  Long,  doing 
different  experiments  have  found  varying  relationships  between  W*  and  E*  (See 
Kantha,  1975). 

The  question  of  why  W*  would  tend  to  be  linearly  related  to  E*  in  many 
but  not  all  cases  may  be  answered  by  closely  examining  the  turbulent  kinetic 
energy  budget  in  the  entrainment  zone. 
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4.4  Turbulent  Kinetic  Energy  Budget  at  the  Density  Interface  and 
the  Development  of  a  Theoretical  Equation  for  Turbulent 
Entrainment  in  the  Presence  of  Mean  Shear 

At  (z  =  -  h),  at  the  top  of  the  entrainment  zone  and  within  the  fully- 
turbulent  mixing  region,  the  turbulent  kinetic  energy  equation  is 


0  =  - 


—  3U  ,  —  9V 
uw  —  +  vw  — 

U   9Z       3Z 


+  bw 


-h 


-h 


dz 


W(P 


P..  +  v 


-  e 


-h 


-h 


(4.9) 


At  (z  =  -  h),  the  turbulent  fluxes  are 


cw 


-h 


Ar  dh 
AC  dt 


(4.10) 


and 


bw 


-h  =  "  *B  oT 


(4.11) 


The  convergence  of  flux  of  turbulent  energy  at  the  interface  is  re- 
sponsible for  the  entrainment  buoyancy  flux,  bw(-h). 


bw"  (-h)  -  I- 


3_ 
3Z 


w(£-  +  I) 

0 


Pn    2 


-h 


The  problem  is  to  estimate  the  time  scale  (re)  required  to  transport  some  of 
the  turbulent  energy  <E>  to  the  vicinity  of  the  entraining  interface. 


9Z 


<r + 1' 


<E> 


-h 


The  mixed  layer  depth  (h)  or  a  length  scale  proportional  to  (h)  is  the 
distance  over  which  turbulent  energy  must  be  transported  by  the  vertical 
component  of  turbulent  velocity  (w).  Therefore,  (xe)  is  taken  to  he  pro- 
portional to  (h)  divided  by  the  rms  vertical  velocity  scale,  /  <w"z>,  giving 
the  entrainment  hypothesis,  equation  (4.12).  This  equation  is  a  refinement 
of  what  Tennekes  (1973)  has  suggested.  The  difference  is  the  use  of 
/  <w"z>  <E>  rather  than  <E>3/2. 
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9_ 
3Z 


w 


(I-  +  I) 


=  mi 


/ 


<w2><E> 


J-h 


(4.12) 


Two  possibilities  are  suggested  for  how  the  mean  shear  at  the  interface 
should  adjust.  The  first  is  according  to  a  gradient-diffusion  model  having 
an  eddy  viscosity  scaling  with  the  integral  scales  of  the  mixed-layer  turbu- 
lence. The  mean  shear  becomes 


3Z 


-  m. 


cw 


i-/<wT> 


(4.13) 


Then  if  local  dissipation  is  negligible  compared  with  the  flux  divergence  at 
the  interface,  equation  (4.9),  using  (4.10)  -  (4.13),  becomes 


m5     lAC|2      au    2 
(— ) 

^dtj 


h  /<wr> 


■  djr 
'df 


m^ 


y<wz>  <e> 


ab  m  +  — u =  o 


(4.14) 


For  those  cases  where  <wz>  is  proportional  to  <E>,  (4.14)  can  be  reduced  to 

simply 


*      p-i       *    mh      *    * 

(W  )2  -  —  -  (W  )  +  —  E  Ri  =  0 

mg         mg 


(4.15) 


in  non-dimensional  form.  Solving  (4.15)  for  W*  gives 


*   Ri 
W  =  — 


2m  k 


Iff 


-m5 


mC  *  * 
—  E  Ri 

mg 


or 


w*  =  «i 


*   r- 


2m 


i  _  /m^* 


5   L 


(4.16) 


where 


4>*  -  mC  ™5  rT*  • 


(4.17) 
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If  <j>*  <  k,   the  binomial  expansion  gives 


/l  -  4  cf>*  =  1  -  2  4>*  -  2  <f>*2  -4  <}>*3  - 


and  therefore  (4.16)  may  also  be  written  as 


W*  =  mC  E*  (1  +  <(.*  +     2  <f>*2  +  .    .    .    .) 


(4.16a) 


For  small  <)>*,  (4.16a)  rapidly  converges  and  only  the  first  term  may  be 
needed.  For  more  general  applications,  where  <w^>  is  not  always  proportional 
to  <E>,  the  equivalent  solution  to  (4.14)  is 


dh 
dt 


m^  /<w2>  <E> 

Rab 


mi+m5<E>|AC|2 

1  +  + 

(hAB)2 


(4.18) 


A  second  possible  way  for  the  mean  shear  at  the  interface  to  adjust  is 
so  as  to  maintain  a  gradient  Richardson  number  of  critical  size,  Ri5.  Then 
in  place  of  (4.13),  the  thickness  of  the  interface  (6)  adjusts  so  that 


3C 
dz 


AB    1 


-h  "  TaC  Ri 


(4.19) 


where 


Rl-  -  ,6AB., 


(4.20) 


In  this  case,  equation  (4.9)  becomes 


M   _L_  \   ar  dh  +  mk  l^l     <E>  -  n 

(1  _  RT~  }  AB  dt  +   h^ ° 

0 


(4.21) 
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or 


dh 
dt 


m^  /<wT>  <E> 
FaB 


(4.21a) 


If  Ri^  is  a  constant  of  order  one,  then  (4.20)  predicts,  independent  of  the 

entrainment  equation  (4.21a), 


6   _1 
h  ~~  Ri 


(4.22) 


4.5  Comparison  with  Moore  and  Long  (1971)  Experiment 

Figure  10  is  a  fit  equation  (4.16)  to  the  results  of  Moore  and  Long 
(1971).  The  constants  (m£)  and  (1TI5)  are  determined  by  the  best  fit.  Notice 
that  in  this  experiment  (E*)"1  happens  to  be  proportional  to  Ri*.  This  is 
probably  due  to  the  experimental  method  and  has  no  general  significance  for 
the  ocean  mixed  layer.  Also,  W*  can  be  interpreted  only  as  a  dimensionless 
entrainment  flux  for  this  experiment  because  the  experiment  involves  two 
turbulent  mixed  layers,  each  entraining  the  other  equally,  resulting  in  a 
motionless  interface,  giving 


VIquah  10.     lk.10n.QXLQ.0Ji  o,ntA.ain- 
mo.nt  cujivq.  v6 .   data  ofi  Uoon.0. 
and  Long    [1971).     An  -inAtab-LLittj 
u>  pn.QdX.oXQd  at  E*  =  0.036.     ThQ 
btopQ  appn.oao.hoj>    (-7)  <X6  E*  do.- 
on.Qai>  qj>  . 


/(39.4E*) 
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w*  = 


-bw(-h)  _  -uw(-h) 
AB  /<£>       AU  /<E> 


The  lowest  order  term  of  (4.18)  agrees  with  the  general  theory  of  Kraus 
and  Turner  (1967)  and  the  experiment  of  Kato  and  Phillips  (1969)  for  those 
instances  where 


/<viT>     <E>  a,  m  iu3  -  u+b+h 


*u*' 


The  higher  order  terms  of  (4.18)  contribute  to  the  entrainment  process 
by  a  feedback  mechanism  that  seems  to  explain  the  large-scale  instability 
observed  by  Moore  and  Long  for  small  Ri*.  Equation  (4.18)  predicts  an 
instability  when  Ri*  <   (Ri*)cr 


where  (Ri*)Cr  =  4  m^5   HaF  '  (4*23) 


For  the  Moore-Long  experiment,  (E*/Ri*)cr  =  0.105,  or 


<Ri*>cr  ■  9-55  m  •  <4-23a» 


Although  this  instability  seems  to  be  possible  in  the  laboratory,  it  is  not 
clear  whether  it  can  ever  occur  in  the  oceanic  mixed  layer  because  the  ex- 
pected turbulent  entrainment  attributable  to  the  turbulent  flux  convergence 
(4.12)  alone  is  sufficient  to  increase  (hAB)  at  a  more  rapid  rate  than  Ri* 
may  decrease. 

The  instability  observed  tends  to  verify  the  gradient-diffusion  hypo- 
thesis, equation  (4.13),  at  least  for  small  Ri*.  However,  the  constant  -Rig 
hypothesis,  equation  (4.19)  seems  to  be  a  better  fit  for  9U/3z(-h)  over  most 
of  the  range  of  E*.  Figure  11  is  a  plot  of  the  observed  (6/h)  versus  the 
values  predicted  by  equation  (4.22). 

Of  course,  the  side  walls  in  flume  experiments  of  this  type  prohibit 
pushing  the  results  too  hard,  but  it  is  fairly  clear  that  the  Moore  and  Long 
experiment  as  well  as  the  others  mentioned  verify  the  lowest  order  term  of 
(4.18a)  as  well  as  (4.21a)— no  matter  how  3C/3z(-h)  is  modeled.  Therefore 
this  shall  be  taken  as  the  entrainment  equation: 
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aBaU 


xlO- 


dh    .,   ^    /<W 

dT  "  hAB 


/<w^>     <E> 


(4.24) 


4.6     Completed  Model   for  Shallow  Mixed  Layers 

Equations   (1.24),    (1.25),   (3.5a,b   )  and   (4.24)  constitute  a  closed  set 
of  five  equations  in  the  five  mixed-layer  variables: 


h  ^-  =  -  AC  |jjr  A  -  1f  <C>h  -  cw(0)    . 


(1.24) 


h  dfB>  __  .  AB  dh  A  +  _M 


dt  dt  p„C 


Qdz  -  bw(0)    . 


(1.25) 


o"P 


-h-6 


1  d 


AC  I2  dh 


2mi 


i  dT  (h  ^>]  - m3u3  +    2~ dT A  - m2/<r>  (<^ 


«ni  3/2 

-  2<v\r>) j-  <E> 


(3.5a) 


(h  <W>)   =  -  ${u*b*  +  AB  ^  A)   +  m2   v^ET  (<q*>     -  2  <w^>   )   -  J-  <E>3/2  • 


2  dt 


dt 


,,        mif/<wT>  <E> 
dT  :  hAB 


(3.5b) 
(4.24) 
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In  addition  to  Q  (z,t),  other  specified  (external)  conditions  are  the  surface 
fluxes  u*b*  and  cw(0). 

Starting  with  initial  conditions  on  <C>,  <B>,  <qz>,  <wz>,  and  h,  the 
model  is  solved  as  an  initial-value  problem  in  time. 

Also  to  be  specified  separately  are  the  initial  profiles  of  B  and  C 
below  the  bottom  of  the  mixed  layer.  These  two  variables  contribute  to  AC 
and  AB,  which  are  important  for  a  deepening  mixed  layer: 

AB  =  <B>  -  B  (-h  -6)  ;  (4.25) 

AC  =  <C>  -  C  (-h  -<5)  .  (4.26) 


5.  BEHAVIOR  OF  THE  EQUATIONS 

5.1  Nondimensional  Form  of  the  Turbulent  Energy 
and  Entrainment  Equations 

Using  the  surface  flux  scales  u*  and  b*,  new  dimensionless  variables  are 
defined: 


* 
H     = 

u*b*  h 

2m3   u*3 

D*    -    . 

hAB     dh 

o    3  dt  ' 


(5.1) 


(5.2) 


2/3 
mi  '     "r 

E*..  =  (— )    2£  ;  (5.3) 


mi2/3 


u* 


<wz> 


E*33  -  (i7)    ^.  (5-4) 

H*  is  the  ratio  of  buoyant  damping  (production)  from  surface  heating  (cool- 
ing) to  wind-stress  production.  P*  is  the  ratio  of  energy  lost  by  entrain- 
ment (potential  energy  increase)  to  wind-stress  production. 

Invoking  the  quasi-steady  state  assumption  for  the  turbulent  energy 
budget,  the  entrainment  and  turbulent  energy  equations  become 
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and 


P*  =  r-E*ii   /^ '  <5-5) 


0  =  l+^r-p2(E*..  -3E33)  /E^T-fCE*^.)3/2  ,        (5.6) 
0  =  -  H*  -  P*  +  n2  (E*..  -  3E*33)  /E*77  -  I  (E^.)37*  ,      (5.7) 


m^         m2 
where     Pi  =  —  and  p2  =  —  .  (5. 8), (5. 9) 


5.2  Determination  of  the  Constants 

The  ratio  m2/mi  =  P2  is  equivalent  to  (18  £j/A)  where  Ui/A)  is  the 
redistribution  to  dissipation  length  scale  ratio  of  Mel  lor  and  Herring 
(1973).  From  boundary  layer  data,  they  suggest 

—  =  0.05  ±  0.01  . 

A 

Hence  p2  is  of  order  one  and  will  be  taken  to  be  equal  to  one  in  this  anal- 
ysis 

I8l1 
P2  =  ~r-^   1  •  (5.10) 


The  ratio  m^/mi  =  Pi  may  be  determined  from  the  asymptotic  case  of  pure 
convection,  H*  ■»■  -  °°  and  Ri*  =  ».  The  equations  yield 


Pi  =  (iTf)  /,E*33/E*1i  t5-11) 


E*. 


33    1         ? 

where  -^ —  =  t  (1  +  irr~   )   for  large  Ri*.  (5.12) 

t  ..   j     jp2 


The  fraction  (r)  of  turbulent  energy  converted  to  potential  energy  by  entrap- 
ment is  found  to  be  0.036  ±  0.031  by  Farmer  (1975)  from  measurements  under 
the  ice.  Using  p2  =  1  and  r  =  0.036,  equation  (5.11)  gives  (pj  =  0.1). 
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P*  may  now  be  solved  as  a  function  of  H*  and  Ri*.  Figure  12  shows  P* 
(H*)  for  Ri*  =  «  and  Ri*  =  1. 

The  one  remaining  constant  needed  to  complete  the  model  is  m3.  This  may 
be  determined  from  the  Kato  and  Phillips  results,  Figure  7,  together  with 
P*(H*  =  0)  from  Figure  12. 


P*(0)  =  hAB   j}JL  (b*  =  0)  =  0.0227 


2m^u 


dt 


(5.13) 


3U* 


hAB  jj£  (b*  ■  0)  =  2.5  u*3  . 


(5.14) 


This  gives  m3  =  55. 

However,  this  value  seems  to  be  quite  high  for  an  oceanic  mixed  layer. 
The  conceptual  model  for  the  mixed  layer  developed  earlier,  equation  (3.3), 
implies  that 


m3  =  |5C0|/u, 


F-iguAe  12.     Ikod^Z  notation  ^on. 
AkaJtlow  mixzd  Jtcuf2SU>:     pot^ntlaZ 
moAQtj  note,  ofa  incA.e.<ue.  dan  to 
2.nt/uiinmzyit  vzamjla  lay&i  Atabit- 
ity. 
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where  5Cq  is  the  excess  surface  and  wind-drive  current 

5C0  =  C  (Z  =  0)  -  <C>  . 


A  more  moderate  value  of  m3  ^  10  is  expected.   In  the  laboratory  experiment, 
the  production  was  enhanced  by  the  side-wall  boundary  layers,  resulting  in  a 


too-large  value  for  m3 


lother  method  of  determining  (m3)  which  is  suggested  by  Fig- 
than  using  neutral  (H*  =  0)  stability  situations  (storm 
)nly  approximately  neutral  by  virtue  of  a  large  u*3  compared 


There  is  another 
ure  12.  Rather 

Situations     are    OFi  i  j     uppiuAimai-ciji     ncunai      uy     vii    luc    ui      a      laiyc     u^        L-UMipc 

with  u*b*h  ),  cases  of  limiting  strong  stability  where  P*  =  0  and  h  =  h, 
could  be  examined 


H*    =  0.4  =  y*^*!i  ,  so 

max       o    a 
2m3u*3 


m3  =  u*b*h    .  (5.15) 

0.8  u  3 


The  effective  buoyancy  flux  u^b*,  given  by  equation  (4.2)  should  properly  be 
determined  by  measuring  Q(Z),  and  the  mixed  layer  depth  h  =  hr  should  be 
ascertained  to  be  in  the  "retreat  mode":  dh/dt  5  0.  This  method  remains  to 
be  attempted. 


5.3  Comparison  With  Earlier  Models 

The  dotted  straight  line  in  Figure  12  is  equivalent  to  the  solution  for 
P*  (H*)  if  equation  (4.1a)  is  calibrated  to  behave  properly  for  both  the  re- 
treating (H*  =  H*max)  and  neutral  (H*  =  0)  situations.  The  greatest  differ- 
ence between  the  two  models  then  occurs  as  H*  becomes  increasingly  negative. 
This  inability  to  calibrate  the  earlier  model  to  both  stable  summer  and 
convective  fall  situations  has  been  demonstrated  by  Thompson  (1974). 

Nearly  equivalent  to  the  Kraus-Turner  type  model  for  entrainment, 
equation  (4.1a),  is  Tennekes's  model  which  suggests  in  effect  that 


Pi  q/o 

P*  =  r  (E*..)3/2  (5.16) 
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rather  than  (5.5).  In  such  a  model,  the  total  turbulent  kinetic  energy  equa- 
tion, the  sum  of  (5.6)  and  (5.7)  is  used: 


0  -  1  -  H*  -  (1  -  j^)  P*  -  (E^-*)372  •  (5.17) 


The  redistribution  term  does  not  appear  at  all  in  (5.17)  but  this  very  term 
is  the  source  of  an  ambiguity  as  H*  +   H*max  =  1.  As  H*  gets  larger,  E*-ji  and 
hence  P*  go  to  zero  for  the  set  of  equations  (5.16)  and  (5.17).  If  rather 
than  using  (5.17),  we  consider  the  component  equations  together  with  (5.16), 
no  unique  solution  exists  for  H*  when  P*  =  0.  Considering  the  budget  for  the 
vertical  component  only  reveals  that  (H*max)  must  be  less  than  one  and  that 
E*-jj  must  be  greater  than  zero  when  P*  =  0.  This,  however,  violates  (5.16). 
The  alternative  equation  (5.5)  does  not  have  this  problem  and  it  suggests  a 
more  consistent  mechanism  for  layer  retreat. 

As  H*  approaches  its  maximum  value,  more  and  more  of  the  vertical  com- 
ponent of  turbulence  is  damped.  Notwithstanding  the  redistribution  effect, 
eventually  redistribution  cannot  keep  up  with  the  combined  losses  due 
to  buoyant  damping  and  viscous  dissipation.  Virtually  all  of  the  energetic 
large-scale  turbulent  motion  is  two-dimensional,  and  the  vertical  flux  of 
turbulent  energy  (-wT)  into  the  interface  ceases  and  therefore  entrainment 
stops.  If  (b*)  is  further  increased  (or  u*  decreased),  (h)  must  retreat 
since  H*  cannot  exceed  H*max  %  0.4.  Locally,  at  the  base  of  the  mixed  layer, 
this  would  correspond  to  the  flux  Richardson  number  (Rf)  surpassing  a  criti- 
cal magnitude  (-  0.2). 

Summing  up,  a  new  mixed-layer  bulk  model  has  been  proposed  that  suggests 
that  the  energy  redistribution  between  vertical  and  horizontal  components  of 
turbulence  must  be  accounted  for  in  a  general  model  used  to  treat  the  entire 
range  of  stability—from  pure  convection  to  layer  retreat.  The  proportion  of 
vertical  turbulent  energy  to  total  turbulent  kinetic  energy,  E*33/E*ii  of 
Figure  13.  varies  most  significantly  in  the  typical  oceanic  range  of  possi- 
bilities, H*  =  0.4  to  *  -  1. 

Layer  retreat  is  predicted  to  occur  for  a  smaller  value  of  stability 
than  is  required  with  the  earlier  models  because  viscous  dissipation  remains 
important  as  long  as  u*  >  0. 

The  model,  as  it  stands,  is  incomplete.  The  above  applies  only  to  the 
relatively  shallow  boundary  layer  associated  with  summer  conditions.  An 
additional  feature  is  required  for  year-round  application  of  the  model 
because  of  the  importance  of  the  rotational  time  scale  in  the  physics  of  deep 
mixed  layers. 
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6.  DEEP  MIXED  LAYERS:  LIMIT  TO  MAXIMUM  DEPTH 

6.1  Limiting  Dissipation  Time  Scale 

The  need  for  a  constraint  upon  the  maximum  depth  of  the  mixed  layer  has 
been  demonstrated.  Niiler  has  pointed  out  that  models  of  the  Kraus-Turner 
type  will  continue  to  deepen  incessantly  during  each  winter  cooling  period. 
This  is  the  unavoidable  consequence  of  requiring  some  fraction  of  the  net 
production  of  turbulent  energy  to  be  allotted  to  increasing  the  potential 
energy  by  means  of  entrainment. 

Thompson  (1974)  also  found  this  to  be  true.   In  doing  a  year-round 
integration,  the  model  required  an  unrealistic  value  for  mean  annual  surface 
heat  flux  yielding  a  net  heat  storage  in  order  to  achieve  a  realistic  cycle. 

Outside  of  the  tropical  latitudes  and  away  from  the  continental  shelves, 
lateral  advection  and  upwelling  are  insufficient  to  balance  this  heat  stor- 
age. The  most  noticeable  failure  of  the  model  is  at  mid-latitude  regions 
having  a  permanent  pycnocline. 

The  model  for  dissipation  is  still 


e  ^ 


<E> 


However,  the  dissipation  time  scale  (te)  was  previously  (in  equation  3.1) 
governed  by  the  depth  of  the  mixed  layer  and  the  turbulent  velocity  scale, 
giving  a  strictly  convective  time  scale: 


Ul>         /<E> 
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As  planetary  rotation  begins  to  rotate  the  mean  shear  profile,  and  hence 
the  geometrical  aspects  of  the  integral  scale,  there  arises  a  new  length 
scale,  associated  with  the  coriolis  time  scale  (f_1). 


l2   =  /<E>  f  . 

It  is  now  increasingly  clear  from  such  studies  as  those  by  Arya  and 
Wyngaard  (1975)  and  Sundararajan  (1975)  that  the  coriolis  time  scale  (f-1) 
plays  an  important  role  in  the  internal  structure  of  the  convective  planetary 
boundary  layer  or  mixed  layer.  The  concern  here  is  more  with  the  bulk  pro- 
perties of  the  region  and  less  with  the  details  within  the  mixed  layer. 
However,  it  is  suggested  that  this  time  scale  has  an  important  role  in  the 
overall  turbulent  energy  budget  above  and  beyond  the  turning  and  decay  with 
depth  of  the  mean  velocity  in  the  classic  Ekman  spiral.  In  reality  the  two 
effects  are  inseparable  because  of  the  link  through  local  shear  production. 

Rather  than  to  simply  replace  the  convective  scale  (&i)  by  the  rota- 
tional scale  {12)9   the  suggestion  is  to  use  £-1  =  i^1  +  i^   ,  giving  equation 
(6.1), 


T 

e 


»-^+f. 


(6.1) 


Replacing  (3.3)  in  the  bulk  (vertically  integrated)  formulation  is  equation 
(6.2) 


i  {E*u)*/2+   Ro"1  (E*.-)  (6.2) 


2m3  u*3 


where 


Ro-*  -  (^)1/3  *  (6.3) 

m3    u* 


and,  as  before, 


E*..  =  (—)2/3   <I> 


11   xm 


3      „2 
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6.2  Nondimensional  Solution  to  the  Entrainment  Function 

Equations  (5.6)  and  (5.7)  will  now  change  to  reflect  the  new  term,  giv- 
ing (6.4)  and  (6.5). 


0  =  1  +  £f*  -  P2  (E*t1  -  3  E*33)  /E*77  -  |  E*..  (/E*77  +  Ro"1).   (6.4) 


0  =   _  H*  -  P*  +  p2   (E*..   -  3  E*33)   /E*77  -  i-  E*..    (/DTT  +  R0-i).    (6.5) 
The  entrainment  equation,    (5.5),  remains  unchanged. 


Pi 


P*  =  -r  E*ii   /^  *  (5-5> 


The  solution  to  P  is  now  a  function  of  two  variables,  and  is  shown  in 


P*  =  P*  (H*,  Ro"1) 


(Figure  14) 


Specifically,  several  aspects  concerning  this  nondimensional  solution 
should  be  observed: 

a.  The  solution  for  Ro-1  =  0  is  identical  with  the  earlier  solution 
(Figure  12).  This  limit  applies  to  both  the  case  of  shallow  (h«u*/f)  mixed 
layer  and  the  case  of  \fery   low  latitude,  <\>  «   u*/|ft|h. 

b.  For  a  retreating  (P*  =  0)  mixed  layer,  hr  is  now  a  function  of  both 
the  Rossby  number  and  the  stability  H*.  Actually,  hr  is  uniquely  determined 
by  the  product  of  these  two  parameters,  making  it  a  function  of 

*    b 
B  =  -2*. 


fu, 


A  fairly  good  analytical  approximation  to  this  function  is  the  hyperbolic 
relationship  (6.6): 
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Entrapment 
Buoyancy  Flux 

-bw(-h)h    _  p* 
3 


Layer 
Stability 

-bw(Q)h      u* 


(Ro)-'=^ 


TIqojkl  14.     GmvwJL  solution  to  zwtAalwmnvvt  and  tuAbulnvvt  kin&tLc  zneAgy 

<iquuoutiovii> . 


(0.5  +  H   )      (1.725  +  Ro"1)  =  1.552 


(6.6) 


which  can  be  put  in  terms  of  hr,  giving  equation  (6.7) 


h  =  1-725  d-  +  0.5L 
nr  -        2 


1  + 


2.76 


1.725  d'  +  5L 


(6.7) 
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where  d '  is  proportional  to  the  Ekman  depth  of  frictional  resistance  and  L' 
is  proportional  to  the  Obukhov  length  scale. 

The  approximate  solution  for  retreat,  equation  (6.7),  has  limiting 
cases  of  d'  ■*  «  and  L'  -*  «  that  are  in  agreement  with  the  earlier  works  of 
Kitaigorodsky  (1960)  and  Rossby  and  Montgomery  (1935). 

L*  =  hp/H*  .  (6.8) 

d*  =  hr/Ro_1  .  (6.9) 

Physically,  the  reason  for  the  cessation  of  entrainment  is  ultimately 
the  same  for  the  convective  boundary  layer  whether  or  not  there  is  signifi- 
cant surface  buoyancy  flux.  Dissipation  that  is  approximately  equally  divid- 
ed among  the  three  components  of  the  turbulent  kinetic  energy  (IF  +  v2  +  w2") 
causes  <w7>to  be  consumed  more  rapidly  than  it  can  be  replaced  by  the  pres- 
sure redistribution  or  "return  to  isotropy"  effect.  When  this  happens,  the 
vertical  flux  of  all  properties,  including  the  turbulence  itself,  ceases  at 
the  depth  h  =  hr.  Buoyant  damping  or  production  modifies  the  overall  budget, 
changing  the  relative  value  of  hrf/u*,  but  dissipation  plays  the  ultimate 
role. 


Buoyant  production  may  be  sufficiently  strong  (H*  <  -0.5)  that  dis- 
si pat ion  cannot  limit  maximum  mixed  layer  depth.  Then  the  production  of 
<w2">  is  too  great  to  be  overcome  by  dissipation  with  increasing  depth. 


c.  A  cyclical  steady  state  is  now  possible.  The  surface  buoyancy  flux 
may  go  through  an  annual  cycle  having  no  net  heat  storage  and  (h)  will  not 
increase  without  bounds  at  the  end  of  the  cooling  phase. 

The  simplest  steady  state  is  the  neutral  stability  point  represented  by 
the  coordinates  (Ro"1,  H*,  P*)  =  (1.38,0. ,0. ) .  This  is  the  situation  for  a 
neutral  boundary  layer,  (bw(0)  =  0),  where  hr  =  1.38  d'. 

A  simple  cycle  having  a  sinusoidal  variation  in  bw(0)  and  a  constant  u* 
is  qualitatively  depicted  by  Figure  15.  This  depicts  the  locus_of  coordin- 
ates for  (Ro-1,  H*)  a_s_  the  cycle  progresses  through  a  heating  (bw(0)  <  0) 
and  then  a  cooling  (bw(0)  >  0)  phase.  The  initial  [0]  and  [4]  points  coin- 
cide at  the  neutral  steady-state  point.  Between  [0]  and  [1]  buoyancy  flux  is 
directed  downwards  and  is  increasing,  and  the  mixed  layer  retreats  to  its 
minimum  depth  which  occurs  at  [1].  Between  [1]  and  [2]  there  is  active 
entrainment  but  the  rate  of  entrainment  is  slow  because  of  significant  buoy- 
ant damping  in  the  stable  mixed  layer.  This  would  correspond  to  the  summer 
warming  period,  culminating  in  maximum  surface  temperatures  at  about  [2]. 
Cooling  begins  at  [2]  and  reaches  a  maximum  at  [3].  Not  until  the  end  of  the 
period  does  the  rate  of  entrainment  become  significant  because  the  heating 
period  [0]-  [2]  causes  a  potential  energy  deficit  in  the  system  which  must  be 
replaced  before  deep-ending  may  proceed  rapidly.  In  the  time  period  just 
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FiguAe.  7  5.     Cyclical.  Atzady  AtaXe.  In  tkz   (H*,Ro   l)   plant. 
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after  [3]  erosion  of  the  thermocline  is  enhanced  by  both  increased  system 
potential  energy  and  by  buoyant  production  which  remains  strong  even  though 
(bw)  is  decreasing. 

Of  course,  the  real  oceanic  system  is  not  in  an  exactly  repeating  annual 
cycle.  However,  year-to-year  variations  resulting  in  long-term  changes  in 
heat  storage  can  be  predicted  accurately  only  by  a  model  that  can  simulate 
the  hypothetical  cyclical  steady-state.  Any  net  annual  buoyancy  flux  across 
the  surface  that  is  balanced  by  advection  is  not  of  concern  in  studying  the 
model  in  the  one-dimensional  form.  In  any  case,  a  large  seasonal  variation 
in  bw(0)--and  u*2--and  the  resultant  cyclical  response  are  dependent  upon  the 
proper  model  response  for  the  more  simple  situation. 


6.3  Simple  Hypothetical  Cases  Demonstrating  the 
Behavior  of  the  Solution 

6.3.1  Shallow  or  f  =  0  (Ro  »  1). 

(1)  Initial  linear  stratification 

(a)  bw  =  0,  u2.  =  constant 

The  result  is  the  same  as  that  predicted  by  the  Kraus-Turner  type  model 
as  well  as  by  the  Kato  and  Phillips  experiment.  See  Figures  16  and  17. 

dji   u^ 

dt  "  hAB 

giving 


t 


1/3 


The  reason  these  models  are  in  agreement  is  that  with  b*  =  0,  u*  is  the  only 
turbulent  velocity  scale. 

<vF->   <*  <E>  a  u£  . 


(b)  Free  convection:  bw(0)  =  constant  >  0,  u|  =  0. 

This  situation  is  out  of  range  for  Figure  13  because  H*  -*  -  °°  .  How- 
ever, if  the  model  equations  are  nondimensional ized  on  the  buoyancy  flux 
velocity  scale  (u^b^h)1/3,  rather  than  on  (u*),  the  problem  may  be  solved. 
In  this  case,  h  =  t*/2,  which  is  also  in  agreement  with  the  Kraus-Turner  type 

models. 
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TIqulAc  7  7.     TcmpeAatuAc  pnofaJLc  at  12- 
hou/i  inteAvalA   fcon.  favc  day*,  6ta/vt- 
ing  uilth  an  inttiat  tinca/i  6  t/icuti  fa- 
cation.     Thcfic  u>  no  iu/ifiacc  heat 
faux  and  a  constant  uxind  i>txa>i>  Lt> 
imposed  at  time  Z2A0. 


:: 


Deg  C 


(c) 


Cyclical  surface  buoyancy  flux: 
constant. 


bw(0)  «  sin  wt,  u*  = 


Figures  18  and  19  show  the  response  of  the  mixed  layer  to  a  simple 
hypothetical  diurnal  heating  and  cooling  cycle.  The  magnitude  of  the  heat 
flux  was  chosen  to  be  particularly  large  in  order  to  accentuate  the  features 
in  a  relatively  short  time  span.  The  initial  stratification  is  the  same  as 
for  the  earlier  neutral  (bw(0)  =  0)run,  l°C/20m.  The  wind  stress  is  also  the 
same,  0.2  dynes. 

Over  the  five-day  period  there  is  an  overall  cooling  trend  because  of 
the  entrainment  process  and  lack  of  net  surface  heat  flux.  However,  the 
daily  maximum  surface  temperature  is  greater  than  the  initial  surface  tem- 
perature for  all  five  days. 

Starting  with  the  fourth  day,  there  is  larger  retreat  corresponding  to 
the  decrease  in  the  Obukhov  length  scale  relative  to  the  value  of  h.  In 
spite  of  midday  retreat,  the  mixed  layer  more  than  recovers  the  temporarily 
decreased  vertical  extent.  Without  eventual  enhancement  of  dissipation, 
deepening  would  continue  without  bound. 

A  comparison  between  Figures  18  and  16  demonstrates  another  important 
difference,  in  addition  to  the  presence  of  a  strong  diurnal  layer  temperature 
change.  Although  there  is  no  net  surface  buoyancy  flux  for  either  case,  the 
situation  having  the  diurnal  cycle  has  a  distinctly  slower  average  rate  of 
deepening.  This  holds  even  for  the  first  three  days  when  there  is  no  re- 
treat. 

The  explanation  for  this  phenomenon  lies  in  the  nonlinearity  of  P*  (H*) 

(Figure  12) . 
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FIquac  18.     Solution  to  hit]    {on  coac  with  cyclical  iunfiacc 
heat  {lux-,   ew[0)   =  0.015  bin   (wt).     CompaAc  with  VIqvjkl  16. 

(d)     General:     u*2  f  0,  bw(0)  f  0 

Since  the  earlier  model  projects  a  linear  P*(H*),  only  two  points  (at 
two  values  of  H*)  may  be  made  to  agree  with  this  model.  For  example,  if 
summer  "retreat"  and  the  neutral  (bw(0)  =  0)  cases  are  coincident  in  both 
models,  then  the  winter  deepening  rate  will  be  greater  for  the  earlier 
model.  This  is  the  mathematical  difference  in  the  basic  responses. 
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VlquJin   7  9.      TrnpoAcutuAQ.  pn.o{ltz  {on. 
cai><L  with  uiitloJL  tinnan.  i>tn.cuti.{i- 
coutlon,  with  cyctlcai  4uA{ace.  hzcut 
{lux.     CompaAe.  mJjfi  FiguAz  17. 


Deg  C 


The  physical  difference  may  best  be  seen  from  the  Tennekes  (1973)  model 

in  which 


-bw 


] 


-h 


<E> 


3/2 


and 


"r  3/2        a       ,   , 

<E>    <*  m  u£  -   u^.b^h  . 


This  model  is  identical  in  response  to  the  Kraus-Turner  model  and  is  suitable 
for  comparison  here  because  it  identifies  the  turbulent  velocity  scale  /<E>. 
Then  retreat  (dh/dt  <  0)  can  occur  only  when  <E>  ■>  0! 

Tennekes  did  not  intend  for  this  formula  to  be  extrapolated  to  the 
strongly  stable  situation.  Instead,  he  suggested  it  as  an  interpolation 
between  the  neutral  and  free  convection  cases—the  typical  atmospheric 
boundary  layer  range.  Nevertheless,  the  above  is  equivalent  to  recent 
oceanic  applications  of  the  Kraus-Turner  type  model.  The  particular  problem 
with  applying  this  method  to  the  stable,  retreating  mixed  layer  is  that  there 
must  still  be  significant  turbulent  energy  (<E>  >  0)  above  z  =  -  hr.  Conse- 
quently, the  fraction  of  turbulent  energy  going  to  increase  the  potential 
energy  cannot  be  independent  of  stability,  H*.  This  explains  the  curvature 
of  P*(H*)  and  the  slower  rate  of  damping  for  the  hypothetical  case  with  an 
oscillating  surface  buoyancy  flux. 
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(2)  Cyclical  surface  buoyancy  flux 

No  1-D  steady  state  is  achievable  because  for  f  =  0,  h  can  be  limited 
only  by  the  Obukhov  length  scale,  therefore  requiring  a  net  downward 
buoyancy  flux  (or  heating).  Then  bw(0)  could  be  balanced  only  by  advection. 
Gill  and  Niiler  (1973)  have  shown  that  advection  is  insufficient  for  lati- 
tudes greater  than  about  15°. 

6.3.2  Convective  planetary  boundary  layer  with  order  one  Rossby  number 

(1)  Neutral  steady  state 

With  zero  surface  heat  flux,  a  steady  state  is  predicted  (see  Fig.  14  at 
the  limit  h  ->  1.38  u*/f ) .  Starting  from  an  initial  linear  stratification  and 
h  <  u*/f,  steady  state  would  be  approached  only  after  a  relatively  long  time. 
This  is  not  likely  to  be  achieved  in  geophysical  flows,  except  perhaps  by  the 
atmospheric  boundary  layer  during  the  polar  winter.  See  Businger  and  Arya 
(1974). 

(2)  Cyclical  steady  state 

This  situation  has  already  been  qualitatively  described  and  illustrated 
by  Figure  15,  but  the  relative  importance  of  the  rotational  and  surface 
buoyancy  flux  scales  needs  to  be  estimated. 

For  the  sake  of  studying  the  relative  response  of  the  mixed  layer,  all 
cycles  will  be  assumed  to  be  repetitive  (in  steady  state).  A  sinusoidal 
surface  buoyancy  flux  and  a  constant  wind  stress  are  used  to  drive  the 
model . 

T  =  pQ  U*2  . 

bw(0)  =  -  |  bw(0)  |  sin  ait  . 

The  layer  response  is  found  to  be  a  function  of  the  parameter  B*,  the 
ratio  of  the  buoyancy  flux  scale  to  the  rotational  scale. 


B*  =  1^(011  .  {6A0) 

f   u*2 

Figure  20  shows  the  cycle  of  h(t),  nondimensionalized  on  u*/f  as  a  function 
of  B*. 

For  large  B*,  which  is  typical  of  the  annual  cycle  for  temperate  oceanic 
regions,  the  low  minimum  mixed  layer  depth  at  cot  =  tt/2  is  attributable  to  the 
strong  influence  of  the  buoyant  damping  corresponding  to  a  small  positive 
Obukhov  length  scale. 
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flgu/tc  20.     Cyclical  6tcady-6tatc 
mixed- lay  qa  n.ci>pon6c  to  6-inu- 
i>oidal  AuA&acc  buoyancy  &lux 
and  constant  wi,nd  6&L2A6. 
t=p0u*2;  HuJ[0)=-\'EuJ[0)\&4.n[tat) . 


In  addition  to  the  change  in  the  range  of  h  with  B*,  there  is  an 
important  change  in  the  shape  of  the  function  h(t).  For  the  case  of  weak 
heating  and  cooling,  the  variation  in  h  is  almost  in  phase  with  bw(0).  In 
the  H*  -  Ro"1  plane  (Fig.  21),  the  locus  of  the  cycle  would  be  small  and 
elongated.  For  increasingly  larger  B*,  the  heat  and  therefore  buoyancy  is 
stored  during  the  "summer"  at  an  increasingly  shallower  depth  (in  comparison 
with  hmax  ^  u*/f).  This  creates  a  substantial  potential  energy  deficit  which 


ViquAQ,  21.  Cyclical  i>tcady-i>tat<i 
in  the  (H*,Ro_1)  plane  ion.  veiy 
hmall  cyclical  buJi^acc  buoyancy 
\Zjjlx   (8*«7). 
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results  in  a  hysteresis  effect  by  retarding  the  entrainment  process.  This 
phase  lag  in  the  heat  storage  naturally  has  important  implications  for  the 
interaction  with  the  atmosphere  for  all  cyclical  time  scales— from  one  day 
periods  of  climatological  importance. 


to 


(3)  Fluctuations  in  wind  stress:  storms 

Fluctuations  in  the  wind  stress,  particularly  in  the  instance  of  storms, 
are  wery   important  in  the  annual  evolution  of  the  mixed  layer  and  the  heat 
storage  of  the  upper  ocean. 

Running  the  model  for  a  simulated  annual  heating/cooling  cycle  as 
before,  but  adding  periodic  "storm"  pulses  of  wind  stress  gives  qualitatively 
little  apparent  change  in  h(t).  (see  Fig.  22).  However,  the  spikes  in  the 
retreat  side  caused  by  the  storms  greatly  enhance  the  storage  of  heat  at 
depth.  The  midsummer  profile  (Fig.  23)  with  the  step-like  remnants  of 
earlier  storms,  has  an  overall  potential  energy  significantly  greater  than 
would  be  the  case  with  no  storms.  This,  in  turn,  will  tend  to  enhance  the 
eventual  deepening,  countering  to  some  degree  the  hysteresis  of  the  heat 
storage  in  the  annual  cycle.  The  annual  range  of  surface  temperatures  (see 
Fig.  24)  is  less  than  it  would  be  without  storms  because  of  the  increased 
heat  storage  at  depth. 
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6.4  Filtering  Effect  of  the  Storage  of  Turbulent  Kinetic  Energy 

The  hypothetical  solutions  shown  in  Figure  20  are  all  for  periods 
(2ir/w)  much  longer  then  the  integral  time  scale  of  the  turbulence,  governed 
by  the  dissipation  time  scale  (xe).  This  made  possible  the  assumption  of 
quasi-steady  state  so  that  the  turbulent  energy  storage  terms,  9/9t  <uju-j>, 
could  be  neglected  in  the  solutions.  However,  fluctuations  in  the  surface 
boundary  conditions  u*2  and  bw(0)  of  sufficiently  short  period  require  the 
consideration  of  these  terms. 

A  wery   simple  model,  which  may  be  solved  analytically,  provides  the 
answers  to  the  question  concerning  the  importance  of  this  unsteadiness  and 

how  it  might  be  treated. 

In  the  simple  model  only  the  two  largest  terms,  shear  production  and 
dissipation,  are  balanced  by  the  time  rate  of  change  of  the  turbulent  kinetic 
energy. 


#  +  ^=^(l  +  asin.t) 


(6.11) 


The  prescribed  shear  production  is  a  constant  plus  a  sinusoidal ly-varying 
component  (a  <  1).  By  defining  dimensionless  time, 


.*  - 


=  t/- 


(6.12) 


VIqixaq.  24.     SuAfcace.  trnpeAatuAn 
T(2=0,t)    faon,  the.  hypothetical  coie 
Mtth  htonmi,'-     CojiH  A,  modeJL  without 
dluAndl  cycZe.. 
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and  turbulent  energy, 


u/*    = 


h<E> 


(6.13) 


(6.11)   is  transformed  to  (6.14), 


W1 

3x^ 


+  y*  =  1  +  a  sin 


0J*T*, 


(6.14) 


where 


to*    =    T        10. 


(6.15) 


For  (a  =  0)  the  steady-state  solution  is  simply  v*  =  1.  For  finite  (a),  the 
complete  solution,  neglecting  the  initial  transient,  is 


y*  =  1  +i 


j*2  +  1 


(sin  w*t*  -  a)*  cos  w*t*) 


(6.16) 


There  is  a  phase  shift  that  increases  with  (u>*),  but  the  important  aspect  of 
(6.16)  is  the  relative  magnitude  of  the  response  associated  with  the  fluctu- 
ating component.  Figure  25  is  a  plot  of  this  versus  to*. 

For  to*  <  *  0.5  (slow  fluctuation  period  compared  with  xe),  the  situation 
of  quasi-steady  state  is  satisfied.  Depending  upon  the  Rossby  number  and  H*, 
this  corresponds  to  a  minimum  fluctuation  period  of  several  hours  up  to  about 
a  day. 


FiguAo.  25.     ReAponie.  o&  the.  mixzd 
layeA  to  fcluctuationA  in  the. 
AuAfaace.  boundary  condition*: 
lilteAlng  d^zct. 
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For  faster  fluctuation  having  a  period  less  than  about  an  hour,  the 
phase  shift  (lag)  in  the  response  approaches  90°  and,  most  importantly,  the 
amplitude  of  the  response  is  negligible.  In  other  words,  the  high  frequen- 
cies are  filtered  out. 

This  becomes  a  possible  cause  of  error  in  using  observed  winds  to  drive 
a  model  having  an  integration  time  step  smaller  or  comparable  with  re.      If 
the  quasi -steady  state  assumption  is  made  to  facilitate  solution,  but  the 
surface  boundary  conditions  are  not  properly  smoothed  or  filtered,  an  incor- 
rect high-frequency  response  will  not  only  be  present  but  will  bias  the  mean 
trend. 


6.5  Interaction  Between  Forcing  Time  Scales:  Modulation  of  the 
Longer-Period  Trend  in  the  Diurnal -Period 
Heating/Cooling  Cycle 

As  already  shown  in  comparing  Figure  18  with  Figure  16,  if  H*  is  per- 
turbed by  a  diurnal-period  heating  cycle,  then  the  mean  P*  over  the  cycle 
will  be  less  than  P*  of  the  mean  H*. 


P*  (H*)  <  P*  'WJ 


(6.17) 


An  atypical  situation  was  used  in  Figure  18  to  demonstrate  the  effect  of 
(6.17)  during  a  relatively  short  time  period.  Over  a  year,  however,  the 
consequences  can  be  significant.  This  is  shown  by  Figures  26-28— case  "B." 
Case  "B"  is  the  same  as  case  "A,"  Figures  22-24,  except  for  the  superposition 
of  a  diurnal  cycle  on  top  of  the  annual  heating/cooling  cycle. 
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Throughout  the  year,  the  principal  effect  is  the  decrease  in  mixed  layer 
depth.  The  decrease  in  upper  ocean  potential  energy,  higher  surface  temper- 
atures, and  lower  heat  storage  at  depth  are  directly  due  to  the  modulation  of 
the  longer-term  entrainment  rate. 

The  basic  shapes  of  the  two  solutions  to  h(t)  are  similar.  Without  a 
noticeable  difference  in  phase,  there  is  the  possibility  that  the  effect  may 
be  absorbed  in  the  constants  of  calibration.   In  other  words,  a  one-day  or 
longer  time  step  might  be  possible  if  computation  speed  were  of  the  essence 
in  a  numerical  model. 


6.6  Simulation  of  a  Real  Case 


An  experiment  conducted  by  Hal  pern  (1974)  provides  a  unique  record  of 
1-D  mixed  layer  evolution  in  that  wind,  currents,  and  thermistor  chain  meas- 
urements were  made  simultaneously  at  one  location  (47°N,  !28°W)  for  the 
period  of  a  month.  During  this  time  one  "strong"  and  several  lesser  storm 
events  occurred.  Figures  29-33  show  the  model  results. 
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FlguAe.  30.     Composition  beArtzzn  obi>eAve.d  and  predicted  tmpeAj&tuAeA  out 

9.6-meXoA  depth. 


Here  the  31-day  records  of  temperature  at  9.6  and  23.4  meters  are  com- 
pared with  the  model  prediction.  The  actual  wind  speeds  are  used  to  compute 
u*  for  the  model,  but  surface  buoyancy  flux  is  rather  crudely  assumed  to  be 
constant  plus  a  constant-amplitude  diurnal  component.  Nevertheless,  the 
gross  features  of  the  record  are  reproduced: 

(1)  A  weak  storm  on  August  11  deepens  the  layer  from  13  to  18 
meters. 

(2)  A  period  of  low  winds  from  August  13  to  August  19  resulted 
in  a  retreating  layer  and  warming  at  the  9.6-meter  depth. 
During  this  period  the  guessed-at  surface  heat  flux  gave  a 
too-shallow  depth  of  retreat  that  can  be  seen  in  the  August 
16  temperature  profile.  Therefore  the  predicted  temperature 
at  9.6  meters  is  too  low  until  the  layer  deepens  suffi- 
ciently on  August  17. 
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Vi.Qan.2.   37.  Compa/uAon  boJM2.2.n  obAQAvzd  and  pfiedicLted  tompoJuaAuJiQA  cut 

23.4-m&£2A  duptk. 

(3)  The  relatively  strong  storm  of  August  20  deepens  the  layer 
to  over  24  meters  and  the  temperature  at  9.6  meters  drops 
considerably.  At  this  time  the  23.4-meter  thermistor  for 
the  first  time  lies  within  the  mixed  layer  and  registers  the 
higher  temperature. 

(4)  Subsequent  warming  and  then  cooling  after  August  30  is  due 
to  another  storm.  This  last  storm  is  nearly  as  strong  as 
that  of  August  20,  but  its  work  is  not  as  obvious  because 
it  starts  with  a  relatively  deep  layer. 

Figure  32  shows  another  possible  model  output--the  mixed  layer  wind- 
driven  current.  As  shown  by  Figure  125  the  mean  current  at  the  base  of  the 
mixed  layer  has  only  a  small  effect  upon  the  turbulent  energy  budget  (even 
for  Ri*  ^  1).  The  depth  of  the  mixed  layer,  on  the  other  hand,  directly 
influences  the  amplitude  of  C.  Therefore,  any  surface  current  prediction 
requires  knowledge  of  h(t). 
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6.7  Evaluation  of  Model  Output 

This  simulation  of  a  real  case  is  not  a  strict  test  of  the  model.  The 
surface  buoyancy  flux  is  not  known  sufficiently  well  for  this  run  to  be  used 
either  for  calibration  purposes  or  as  a  test  of  model  behavior. 

Inasmuch  as  the  model  was  designed  to  have  general  applicability  over 
all  seasons  (wide  range  of  values  for  H*),  the  proper  calibration  and  testing 
procedure  requires  an  integration  over  a  full  year.  This  remains  to  be  done 
at  a  future  date  under  a  more  operationally-oriented  program  of  research. 

Figure  32  shows  a  measure  of  the  sensible  heat  flux  at  the  surface  (by 
a  bulk  parameterization)  versus  time  for  the  first  half  of  the  total  period 
of  observation.  (Air  temperature  was  not  available  after  the  storm  of  21 
August.)  Because  relative  humidity  was  not  available,  the  relatively  larger 
latent  heat  flux  could  not  be  computed.  Because  the  latent  heat  flux  is  very 
roughly  proportional  to  the  sensible  heat  flux  (assuming  an  approximately 
constant  Bowen  ratio),  Figure  32  indicates  the  large  magnitude  of  the  net 
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diurnal  heat  flux  compared  with  the  mean  daily  values.  The  mean  daily  heat 
flux  and  the  diurnal  component  were  assumed  to  be  constant  over  the  entire 
period.  This  was  obviously  not  true. 


Nevertheless 
does  have  didacti 
sponse  to  realist 
particular  case  i 
ment  and  retreat 
tant  in  determini 
guess  at  bw(t)  wa 
the  diurnal  cycle 
integration  using 


,  this  particular  run,  as  illustrated  by  Figures  29  to  33, 
c  value  in  demonstrating  the  character  of  the  model  re- 
ic  boundary  conditions.  The  overall  conclusion  for  this 
s  that  the  wind  stress  controls  the  timing  of  the  entrain- 
events.  The  degree  of  surface  heating  or  cooling  is  impor- 
ng  the  extent  of  the  retreat  or  the  deepening.  Although  the 
s  only  marginally  satisfactory,  an  important  result  is  that 

prescribed  gives  a  much  better  prediction  than  does  an 

a  constant  heat  flux  alone. 


7.  CONCLUSIONS 

This  model  for  the  ocean  mixed  layer  consists  of  a  number  of  hypotheses 
and  assumptions  regarding  the  mechanical  energy  budget,  originally  shown  in 
Figure  1.  In  arriving  at  these  hypotheses,  mean  turbulent  field  techniques 
were  used,  but  this  was  done  without  losing  the  simplicity  of  the  bulk  model 
concept.  The  particular  details  of  the  structure  of  the  mean  fields  within 
the  boundary  layer  were  not  of  highest  priority  as  a  model  output.  In  any 
case  detailed  geophysical  data  revealing  this  structure  is  unavailable  to 
test  such  a  prediction.  Of  highest  importance  in  this  model  was  the  ability 
to  predict  the  year-round  evolution  of  the  ocean  mixed  layer  depth,  together 
with  the  bulk  properties.  There  are  several  models,  including  those  by 
Wyngaard  et  al.  (1974)  and  Sundararajan  (1975),  which  provide  structural 
predictions  quite  nicely,  but  these  boundary  layer  models  are  not  directly 
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concerned  with  a  changing  layer  height.  An  entrainment  model  is  necessary  to 
provide  boundary  conditions  at  the  moving  density  interface. 

The  atmospheric  and  oceanic  boundary  layers  present  nearly  identical 
problems.  However,  an  important  difference  arises  in  the  surface  heat  flux 
boundary  condition.  The  atmospheric  boundary  layer  is  predominantly  unstable 
because  of  minimal  absorption  of  solar  radiation  within  the  layer.  On  the 
other  hand,  most  of  the  solar  radiation  is  absorbed  near  the  ocean  surface 
and  hence  turbulent  heat  flux  downward  in  the  oceanic  boundary  layer  is  about 
the  same  as  flux  upward  in  the  course  of  an  annual  cycle.  Therefore,  the 
stable  situation  is  of  relatively  greater  importance  in  the  ocean.  As  shown 
by  Figures  12  and  14,  the  nonlinearity  of  P*  (H*)  is  greatest  for  H*  >  0, 
the  stable  region.  The  treatment  of  the  oceanic  case  must  be  made  with 
emphasis  upon  this  point  since  the  linear  extrapolation  of  the  unstable 
situation  will  not  suffice. 

In  addition  to  the  nonlinear  dependence  upon  stability,  a  Rossby  number 
dependence  for  the  entrainment  rate  has  been  found  in  the  nondimensional 
solution,  Figure  14.  This  makes  possible  a  cyclical  steady  state  for  the 
boundary  layer  without  resorting  to  unrealistic  values  of  upwelling  or 
lateral  advection  in  a  long-term  integration. 

In  the  short  term,  on  the  order  of  days,  the  upper  ocean,  even  at  higher 
latitudes,  exhibits  significant  three-dimensional  baroclinic  activity.  A 
case  in  point  is  well  documented  by  Gregg  (1975)  in  which  only  lateral  advec- 
tion can  explain  the  change  in  heat  content  after  the  passage  of  a  storm.  In 
doing  the  ensemble  or  horizontal  averaging,  such  phenomena  as  well  as  the 
short-period  internal  waves  are  treated  as  noise  in  the  essentially  one- 
dimensional  mixed  layer  dominated  by  the  strong  vertical  turbulent  fluxes  of 
heat,  salt,  and  momentum.  This  does  not  eliminate  the  practical  difficulties 
of  analyzing  and  simulating  actual  observations  made  at  a  single  point,  as 
shown  by  the  case  using  real  winds  to  drive  the  model. 

The  importance  of  shorter-period  fluctuations  in  the  surface  buoyancy 
flux  in  modulating  the  longer-term  response  has  demonstrated  the  need  to  know 
the  typical  daily  heating/cooling  cycle:  Solar  radiation,  evaporation,  con- 
duction, and  back-radiation  as  a  function  of  season  and  geographical  coordin- 
ates. The  detailed  features  of  this  diurnal  cycle  were  not  of  importance  in 
conducting  a  qualitative  evaluation,  but  accurate  quantitative  results  from 
an  operational  model  would  be  impossible  without  them.  For  the  same  reason, 
the  radiation  absorption  function  should  be  known.  This  last  suggestion  may 
not  be  possible  without  the  incorporation  of  a  primary  productivity  model  in 
which  case  the  physics  becomes  directly  dependent  upon  the  biology! 

Some  of  the  most  fruitful  applications  of  this  model  lie  in  the  link-up 
with  other  models.  The  region  below  the  mixed  layer  has  been  over  a  short 
time  period  successfully  treated  as  a  zero-flux  "quiescent  abyss."  This  was 
done  primarily  for  simplicity  when  the  point  of  focus  was  the  overlying 
turbulent  boundary  layer.  The  first  natural  joining  of  models  is  between 
such  a  boundary  layer  model  with  an  ocean  general  circulation  model.  But  why 
stop  there?  An  atmospheric  general  circulation  model  may  be  combined  with 
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the  ocean  general  circulation  model --but  only  through  two  boundary  layer  or 
mixed  layer  models.  It  has  been  demonstrated  here  that  the  sea-surface  tem- 
perature, which  is  so  important  for  weather  prediction,  cannot  be  predicted 
with  accuracy  without  the  proper  consideration  of  the  ocean  mixed  layer. 
Climate  modeling  must  also  pay  heed  to  the  oceanic  boundary  layer.  Anomalous 
events  in  the  evolution  of  the  mixed  layer  may  induce  a  delayed  response  by 
the  atmosphere  when  the  mixed  layer  again  deepens  sufficiently. 

The  interaction  between  diurnal  and  annual  time  scale  events  has  been 
particularly  emphasized,  but  other  time  scales  even  longer  than  a  year  pose 
some  interesting  possibilities,  particularly  for  inactive  models. 


8.  SUMMARY  OF  MODELED  EQUATIONS 
ENTRAPMENT  BUOYANCY  FLUX 

2 


m^  <w2>  <E> 
■bw  (-h)  =  £ (1) 


BUDGET  FOR  HORIZ.  COMPONENTS  OF  TURB.  KINETIC  ENERGY 


Id     /.      -y—, — ?■  x                 o       h  bw  (-h)  !r>s 

2jt   (h  <uz  +  vz>)   =  m3  u*3   2  R]»  (2) 


2m! 
-m2   (<E>  -  3  <wz>)     <E>"2  -  -T-  (<E>"2  +  fh)   <E> 


BUDGET  FOR  VERTICAL  COMPONENT  OF  TURB.  KINETIC  ENERGY 


^d_(h  <w7>)  =  llM^hl  +  h_M0i  (3) 


+  m2    (<E>  -   3    <wz>)   <E>2  -  ^-  (<E>2  +  fh)   <E: 


CONSERVATION  OF  MEAN  BUOYANCY  AND  MEAN  MOMENTUM 
h 


o 


^JS>  -  Bw  (-h)  -  Bw  (0)  + -SS-   /   Q  dz         (4) 

dt  pocp  J.h 
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h  ^^  =  cw  (-h)  -  cw  (0)  -  if  <C>  h  (5) 


'JUMP  CONDITIONS"  AT  BOTTOM  OF  MIXED  LAYER 


-  CW  (-h)  =  AC  jj: 


-  bw  (-h)  =  AB  ^ 


NOTATION 


L/2     L/2 
C  =  U  +  i  V  (")  =  I™      L.  j  j        (  )  dxdy 


-L/2     -L/2 


E  =  u2  +  v2  +  w2 


u*2  =  |cw  (0)|  <(  )>  =  £  f   (  )  dz 

-h 
B  =  (p0  -  J)   g/p0 

B  =  P0_1  3P/3T 
hAB 


Ri*  = 


AC 
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